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ABSTRACT 

This  program  is  based  on  the  analytical  work  contained  in  Report 
APCRL  65-579  (Martin  Report  ER  13950,  Ref.  1).  In  addition,  some 
analytical  methods  not  covered  in  the  above  report  are  presented  in 
the  appendix  of  this  report. 

The  operating  modes,  general  features  and  accuracy  of  the  program 
are  discussed.  Operating  instructions  and  input/output  descriptions 
and  definitions  are  provided.  All  symbols  used  in  the  program  are 
listed  and  defined.  Flow  charts,  descriptions  and  explanations  of  the 
program  and  subroutines  are  also  included. 

The  program  is  written  in  Fortran  IV  and  machine  language  (MAP). 
Double  precision  is  used  extensively. 
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FOREWORD 

This  computer  program  is  the  result  of  research  which  was 
performed  for  the  Data  Analysis  Branch  (CRMXA),  Technical  Services 
Division  at  AFCRL,  USAF,  L.  G.  Hanscom  Field,  Bedford, 
Massachusetts.  The  contractor*  s  report  number  is  ER  14226. 
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I.  INTRODUCTION 


Orbit  determination  from  satellite  observations  is,  essentially, 
a  process  in  which  a  theoretical  orbit  is  fitted  through  the  observations 
in  a  manner  which  minimizes  some  error  function.  The  theoretical 
orbit  is  represented  by  the  mathematical  model,  and  the  best  fit,  in 
the  present  program,  is  defined  as  one  which  minimizes  the  variance. 
The  analytical  foundations  of  the  Minimum  Variance  Method  employed 
in  this  program  are  presented  in  Ref.  1,  which  also  includes  the 
details  of  the  mathematical  model  of  the  dynamical  system.  The 
present  report  deals  with  the  computer  program  itself. 

The  program  consists  of  two  main  parts:  orbit  estimation  or 
filtering  routine  and  ephemeris  computation  routine.  The  analytical 
treatment  of  the  filtering  routine  is  reported  in  Ref.  1,  which  includes 
some  of  the  equations  utilized  in  the  ephemeris  computation.  The 
specific  methods  used  in  the  ephemeris  computation  are  presented  in 
the  appendix  of  this  report,  which  also  includes  the  development  of 
additional  equations  used  in  the  filtering  routine,  but  not  reported  in 
Ref.  1.  It  must  be  pointed  out  that  although  the  filtering  and  the 
ephemeris  routines  are  two  separate  subprograms,  they  are  inter¬ 
connected  with  extensive  logic,  which  enables  the  program  to  perform 
both  operations  simultaneously  when  required.  This  assures  greater 
accuracy  of  the  estimated  orbit  and  saves  considerable  computer  time. 

The  accuracy  of  the  estimated  orbit  is  dependent  on  three  main 
factors:  the  observations,  the  filtering  process,  and  the  mathematical 
model  of  the  dynamical  system.  The  accuracy  and  frequency  of  the 
observations  are,  of  course,  outside  the  control  of  the  program  but  are, 
nevertheless,  some  of  the  most  important  factors.  One  of  the  most 
difficult  problems  is  the  screening  of  illegitimate  observations  when 
the  interval  between  observations  is  large  and/or  when  the  total 
number  of  observations  is  small.  This  is  because,  under  the  circum¬ 
stances,  there  is  no  real  basis  for  a  rejection  criterion.  To  cope  with 
such  problems,  a  rather  elaborate  rejection  technique  was  developed 
uhich  is  outlined  in  Ref.  1  and  further  discussed  in  Sections  C  and  D  of 
the  Appendix. 

A  special  program  was  developed  to  test  the  filtering  method.  It 
was  established  that  for  normally  distributed  observation  errors,  the 
filtering  process  converges  to  the  true  values  of  the  orbital  elements 
within  the  numerical  accuracy  of  the  computer  by  processing  about  50 
to  100  observations.  The  convergence,  however,  is  asymptotic,  and 
gains  in  accuracy  are  small  in  the  latter  phase  of  filtering. 

The  third  factor  which  affects  the  accuracy  is  the  mathematical 
model  of  the  dynamical  system.  In  the  present  program,  it  includes 
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Jacchia*  s  1964  model  atmosphere’  and  six  zonal  harmonics.  The 
equations  for  additional  harmonics  and  solar  and  lunar  perturbations 
have  been  developed  in  Ref.  1  and  can  be  easily  included  in  the  present 
program.  However,  it  was  considered  that  a  significant  increase 
in  computing  time  resulting  from  their  inclusion  is  not  justified  under 
the  present  circumstances.  The  higher  order  perturbations,  however, 
are  included  in  special  purpose  programs  (Ref.  3). 


n.  PROGRAM  DESCRIPTION 
A.  MODES  OF  OPERATION 

The  program  is  designed  to  perform  two  main  ftinctions:  Obtain  a 
best  estimate  of  the  orbit  and  compute  the  ephemeris.  Thus,  the 
program  could  be  separated  into  two  parts.  In  actuality,  the  two  parts 
are  interconnected  in  the  sense  that  the  ephemeris  is  computed  con¬ 
currently  with  the  filtering  whenever  the  two  time  periods  coincide 
and  the  filtering  is  performed  in  the  forward  direction  for  the  last 
time. 

The  program  can  be  operated  in  three  modes.  The  first  is  a  single 
filtering  mode.  This  mode  will,  normally,  be  used  when  the  initial 
estimate  of  the  orbit  and  the  system  observation  errors  are  well 
known.  This  can  happen  in  restarting  the  filtering  with  orbital  elements 
obtained  from  a  previous  filtering  and  with  observations  for  a  subsequent 
time  period.  If  the  ephemeris  is  required  from  a  time  previous  to  the 
input  values,  the  program  will  integrate  backward  to  the  first  required 
ephemeris  time,  then  integrate  forward  while  computing  ephemeris 
until  the  first  observation  is  encountered.  From  this  point,  filtering 
and  ephemeris  computation  is  done  concurrently  until  the  last  observa¬ 
tion  is  encountered.  If  the  ephemeris  is  required  past  the  last  obser¬ 
vation,  there  is  an  option  for  smoothing  (see  Section  n-D)  wherefrom 
the  integration  and  ephemeris  computation  proceeds  to  the  final  time. 

It  is  obvious  that  in  order  to  utilize  this  mode  of  operation,  the  initial 
estimate  of  the  orbit  must  be  good  if  the  ephemeris  is  required  for 
periods  where  the  filtering  has  not  improved  the  estimate  to  the  desired 
degree.  The  specification  of  the  time  period  for  the  ephemeris  com¬ 
putation  is  arbitrary  but,  of  course,  must  be  kept  within  reason.  A 
schematic  illustration  of  the  mode  is  shown  in  Fig.  1 . 


It  must  be  pointed  out  that  the  initial  values  of  the  orbital  elements 
can  be  given  for  any  time,  but  for  practical  reasons  must  be  close  to 
the  time  of  the  first  observation  or  the  first  ephemeris  time  to  avoid 
excessive  integration. 
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FIGURE  I.  FIRST  MODE  OF  OPERATION 


The  second  mode  of  operation  involves  triple  filtering  and  is  used 
whenever  the  standard  deviations  of  the  observations  for  the  system 
are  not  known  and/or  the  initial  estimate  of  the  orbit  is  not  known  with 
sufficient  accuracy. 

A  number  of  observation  sets  (up  to  200  cards)  can  be  specified 
which  will  be  filtered  three  times.  This  is  done  in  order  to  estimate 
the  standard  deviations  of  the  observations  for  the  system  and  to 
obtain  a  good  estimate  of  the  orbit  in  case  the  ephemeris  output  is 
required  in  a  region  close  to  the  first  observations.  The  process  is 
started  with  the  initial  estimates  of  the  orbital  elements  and  system 
observation  errors  by  filtering  forward  the  specified  number  of 
observation  cards  (normally  20  to  40  or  one  to  two  days  of  observa¬ 
tions).  At  the  end  of  the  first  filtering,  an  improved  estimate  of  the 
system  standard  deviations  and  the  rejection  criterion  (see  Appendix, 
Sections  C  and  D)are  computed,  and  the  filtering  is  continued  backward. 
The  process  is  repeated  during  and  after  the  backward  filtering.  If 
ephemeris  is  required  prior  to  the  first  observation,  the  orbit  is  inte¬ 
grated  backward  to  the  first  ephemeris  time  and  the  ephemeris  computed 
while  integrating  forward.  The  further  process  is  similar  to  the  corre¬ 
sponding  part  of  the  first  operating  mode.  An  illustration  of  the  mode 
is  shown  in  Fig.  2. 

The  rejection  process  is  not  shown  in  the  diagram,  but  is  done 
whenever  filtering  is  performed.  It  was  established  that  no  noticeable 
gains  in  accuracy  were  achieved  by  filtering  more  than  three  times. 

It  must  be  pointed  out  that  the  rejection  criterion  is  continuously  being 
tightened  during  the  filtering  and  not  Just  at  the  end  of  each  filtering 
process  as  in  the  case  with  the  Least  Squares  Method  in  a  multiple 
filtering  mode. 
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FIGURE  2.  SECOND  MODE  OF  OPERATION 


The  third  mode  of  operation  does  not  involve  filtering  and  is  used 
only  for  predicting  the  ephemeris  (including  observations)  from  given 
initial  conditions.  Thus,  it  can  be  used  in  preflight  orbit  analysis. 
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FIGURE  3.  THIRD  MODE  OF  OPERATION 


As  in  the  first  two  modes,  the  first  and  last  ephemeris  times  are 
arbitrary  in  relation  to  the  time  of  the  initial  values. 
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B.  OBSERVATION  ERRORS 


To  (Ully  utilize  the  capabilities  of  the  Minimum  Variance  Method, 
it  is  necessary  to  know  the  standard  deviations  of  the  observations, 
which  are  used  to  obtain  the  weighting  matrix.  In  defining  the  standard 
deviations,  however,  two  factors  must  be  considered.  First,  the 
mathematical  model  of  the  actual  dynamical  system  is  not  perfect. 
Therefore,  in  fitting  the  theoretical  orbit,  the  errors  in  the  mathemat¬ 
ical  model  will  be  accommodated  as  observation  errors.  Secondly, 
there  are  unknown  bias  errors  in  individual  observation  systems 
(station  locations,  etc. ).  The  theory  of  the  Minimum  Variance  Method 
is  based  on  a  normal  or  Gaussian  error  distribution.  The  determina¬ 
tion  of  the  bias  errors  simultaneously  with  the  orbital  elements, 
although  possible,  is  not  practical  for  a  general  type  orbit  determination 
program,  because  the  error  functions  and  bias  constants  will  have  to 
be  defined  separately  for  each  individual  station.  Therefore,  the 
approach  taken  in  the  present  program  is  based  on  the  standard  devia¬ 
tions  for  the  entire  system  instead  of  the  individual  stations.  Since, 
in  all  probability,  the  observation  biases  (including  locations)  for  all 
observing  stations  will  not  be  biased  in  the  same  direction,  the  system 
bias  will  become  more  negligible  as  the  number  of  stations  increases. 
The  standard  deviations  for  the  system  now  will  include  the  biases  of 
the  individual  stations,  but  they  can  be  considered  random  for  the 
system.  The  system  standard  deviations  (which  include  the  bias  errors) 
can  be  estimated  by  successive  approximations  in  a  multiple  filtering 
process.  This  is  done  in  the  present  program  when  operating  in  the 
second  mode  or  triple  filtering. 

Initally,  an  estimate  is  made  of  the  system  observation  errors. 

This  estimate  is  used  in  the  first  forward  filtering.  In  addition,  a 
number  of  sigmas  are  inputted  which  are  used  for  the  rejection  criterion 
(Appendix,  Section  D)  during  the  first  filtering.  Since  the  orbit  and 
the  system  standard  deviations  are  not  known  accurately,  the  number 
of  sigmas  for  the  first  filtering  is  usually  large  (50  to  100).  Thus 
only  grossly  inaccurate  observations  will  be  rejected  during  the  first 
filtering.  At  the  end  of  the  first  filtering,  better  estimates  of  the 
system  standard  deviations  are  computed  based  on  the  filtering  results. 
Also  the  number  of  sigmas  for  the  rejection  criterion  is  now  computed. 
This  number  will  usually  be  three,  unless  the  initial  estimates  of  the 
system  standard  deviations  have  been  considerably  underestimated, 
in  which  case  the  number  of  sigmas  will  be  larger  to  avoid  rejection 
of  legitimate  observations.  These  estimates  are  used  for  the  second 
or  backward  filtering,  at  the  end  of  which  new  and  improved  estimates 
of  these  quantities  are  obtained,  which  now  are  used  for  the  third  or 
forward  filtering  phase. 

During  the  three  filterings,  the  criterion  for  rejection  of  observa¬ 
tions  is  continuously  tightened  as  the  orbit  becomes  known  with  higher 
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accuracy.  However,  the  improvement  in  accuracy  of  the  orbit  becomes 
small  during  the  latter  pu?ses  of  filtering  as  it  is  strongly  dependent 
on  the  restraining  accuracy  of  the  observations.  Thus,  the  criterion 
for  rejection  is  essentially  constant  during  the  last  filtering  phase. 


C.  ITERATION 


Unfortunately,  in  practice,  there  are  many  cases  when  observations 
are  infrequent  and  may  involve  intervals  of  several  days.  If  the  filtering 
has  not  progressed  at  time  t^  to  the  point  .vhere  the  orbit  is  known 

with  sufficient  accuracy,  the  deviation  of  the  estimated  orbit  from  an 
observation  at  time  t^  +  j  after  a  long  interval  may  be  large.  Assuming 

a  legitimate  error  in  the  observation  at  t^  + 1 ,  tide  may,  normally, 

indicate  a  perfectly  legitimate  error  in  the  estimate  of  the  orbit  at 
time  t^.  Since  the  filtering  method  is  based  on  linear  assumptions. 


this  large  discrepancy  between  the  observation  and  the  estimated  orbit 
will  certainly  violate  the  linearity  and  consequently  affect  the  computed 
corrections  and  thus  degrade  the  accuracy  of  the  estimated  orbit.  On 
the  other  hand,  because  of  the  nature  of  orbits,  the  large  discrepancy 
at  time  t^  ±  t  is  caused  by  a  relatively  small  error  in  the  orbital 


+ 1 

elements  at  time  tk. 
one  sigma  accuracy. 


This  error  at  t^  is. 


normally,  well  within  the 


On  the  basis  of  these  considerations,  a  method  was  developed  to 
cope  with  this  problem.  After  computing  the  corrections  of  the  posi¬ 
tion  and  velocity  vector  for  a  ne^v  observation  at  time  t^  +  j  by  use  of 

the  weighting  matrix,  the  corrections  are  tested  against  a  criterion. 
If  they  exceed  the  allowable  limits,  an  iterative  process  is  begun 
whereby  the  corrections  at  t^  +  ^  are  transferred  to  t^  by  use  of  the 

state  transition  matrix 


x(‘k)  *  ♦  (tfc.  tk  +  i>  X  (tk+1)  (1) 

or  its  inverse 

x(tk)  ■  *’1(tk+r  tk)x<tk+i>  <2> 


whichever  is  more  convenient.  The  corrections  to  the  orbit  are  made 
at  tk  and  the  orbit  integrated  again  to  t^  +  If  the  orbital  process 

would,  indeed,  be  linear,  the  iterated  orbit  at  tk  +  ^  would  correspond 
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to  the  linearly  corrected  orbit.  Since  it  is  not.  the  iterated  orbit 
will  be  different.  The  process  is  convergent  and  is  repeated  until  the 
corrections  at  t^  +  ^  become  tolerable,  which  usually  occurs  in  less 

than  three  iterations.  The  method  has  proved  itself  to  be  very  valuable 
in  difficult  cases. 


D.  SMOOTHING 

The  estimated  orbit  is  not  continuous  as  determined  by  the 
Minimum  Variance  Method  in  the  form  of  discrete  position  and  velocity 
vectors  at  the  observation  times.  The  discontinuities  are  represented 
by  the  corrections.  However,  by  using  the  filtering  methods  outlined 
in  Section  II-A.  the  discontinuities  are.  normally,  negligible  during 
the  ephemeris  computation  phase.  Therefore,  smoothing  for  this 
purpose  is  not  considered  necessary  and  thus  considerable  savings 
in  computing  time  can  be  achieved.  However,  in  predicting  the  orbit 
for  a  future  period  of  time,  every  available  means  should  be  utilized 
to  improve  the  final  position  and  velocity  vector,  since  it  becomes 
the  sole  source  of  information  for  the  future  orbit.  In  the  present 
program,  smoothing  is  employed  for  this  purpose.  The  final  orbit 
is  smoothed  through  the  last  six  points  which  are  at  least  10  min. 
apart  and  for  which  the  total  time  span  does  not  exceed  1.  5  days. 

These  values  were  chosen  for  practical  and  numerical  reasons.  The 
theory  is  developed  in  Section  B  of  the  Appendix.  The  prediction  for 
a  time  previous  to  the  first  observation  is  not  based  on  a  smoothed 
value  and,  therefore,  will  be  somewhat  less  accurate  depending  on  the 
accuracy  of  the  estimate. 


E.  NUMERICAL  METHODS 

In  dealing  with  matrices  of  even  moderate  order  (the  largest 
matrix  encountered  in  the  present  program  is  7x7)  and  continuously 
processing  a  large  number  of  data,  it  is  seldom  that  certain  numerical 
problems  do  not  arise.  The  main  difficulties  in  this  type  of  program 
are,  usually,  connected  with  matrix  inversion  and  the  degradation  of 
the  numerical  values  due  to  the  accumulation  of  round-off  errors. 

The  analytical  and  numerical  methods  used  in  the  present  program 
virtually  eliminate  any  difficulties  in  matrix  inversion  (see  Ref.  1, 
Section  V-D).  This  has  been  borne  out  by  extensive  experience,  which 
has  shown  no  evidence  of  such  problems. 

The  second  problem,  connected  with  the  degradation  of  the  numerical 
values,  however,  is  always  present  in  dealing  with  a  large  amount  of 
data  and  extensive  matrix  operations.  Therefore,  it  has  been  necessary 
to  introduce  certain  numerical  manipulations  to  assure  a  continuous 
operation  of  the  program. 
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One  of  the  problems  concerns  the  covariance  matrix.  Theoretically, 
the  covariance  matrix  should  always  be  symmetrical  and  the  diagonal 
elements  should  always  be  positive,  hi  practice,  however,  the  covari¬ 
ance  matrix  becomes  unsyxnmetrical,  and  occasionally  a  diagonal  element 
may  assume  a  negative  value  as  a  result  of  round-off  error  accumula¬ 
tion  in  extensive  matrix  operations. 

To  cope  with  the  symmetry  problem,  the  program  employs  a  process 
whereby  the  covariance  matrix  is  symmetrized  by  equalizing  one  side  of 
the  diagonal  to  the  other  side  after  every  major  operation.  This  is  faster 
than  averaging  and  the  error  introduced  is  negligible,  since  only  the  last 
one  or  two  places  of  an  element  are  normally  unsymmetrized  in  a  par¬ 
ticular  operation,  as  borne  out  by  experience. 

The  negative  diagonal  elements  can  occur  in  the  process  of  trans¬ 
forming  the  covariance  matrix  over  a  long  time  period.  This  implies 
extensive  correlation  of  the  elements.  The  remedy  in  such  case  is  to 
strip  the  matrix  of  the  covariance  elements  and  transform  only  the 
diagonal  matrix.  This  is  done  in  the  program.  This  case  is  more  likely 
to  occur  at  the  beginning  of  the  filtering  process  when  the  initial  estimate 
of  the  orbit  is  rather  poor. 

Another  case  where  negative  diagonal  elements  may  appear  in  the 
covariance  matrix  is  in  the  updating  process.  In  particular,  this  is  more 
likely  to  occur  in  cases  where  accurate  measurements  are  introduced  at 
a  point  where  the  orbit  is  relatively  poorly  known.  In  this  case,  updating 
involves  the  substruction  of  two  matrices  of  equal  magnitude.  This 
problem  can  be  coped  with  by  assuming  the  previous  updated  value  for 
the  particular  element,  or  by  some  other  reasonable  method. 

It  must  always  be  remembered  that,  in  practice,  the  covariance  matrix 
is  not  an  array  of  absolute  numbers  and,  therefore,  can  be  dealt  with  ac¬ 
cordingly.  The  occasional  occurrence  of  the  negative  diagonal  elements 
in  covariance  matrix  operations  itself  proves  this  point  since  theoretically 
it  should  not  happen. 

The  numerical  problems  discussed  above  can  be  alleviated  by  carrying 
more  digits,  in  other  words,  using  double  precision  in  the  matrix  opera¬ 
tions,  This  is  done  in  the  present  program  wherever  necessarv. 


P.  COMPUTATION  OP  EPHEME  RE3 

The  second  function  of  the  program  is  to  compute  the  required  ephem- 
eris  and  related  quantities.  As  mentioned  before,  the  ephemeris  is  com¬ 
puted  concurrently  with  the  last  forward  filtering,  when  the  two  time 
periods  coincide.  In  addition,  it  is  computed  outside  the  filtering 
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region  whenever  required.  The  initial  and  final  times  are  input  values, 
as  is  the  computation  interval. 

In  addition  to  the  regular  printout  times,  up  to  20  odd  printout  times 
can  be  specified. 

The  quantities  computed  and  printed  out  are:  date  and  time,  revolu¬ 
tion  number,  position  and  velocity,  Greenwich  mean  sidereal  time, 
classical  orbital  elements,  observations  and  their  time  rates  for  a 
specified  number  of  stations.  Six  options  are  available  for  computing 
the  desired  quantities.  Fbr  more  detail,  see  Input  and  Output  sections. 

Output  can  be  specified  as  regular  printout,  binary  tape,  and  binary 
coded  decimal  (BCD)  tape. 


W0* 


t  ***.  *V  r  * 


mmtmm  wmm  i 


ill.  PROGRAM  INPUT 


A.  NOTE  ON  THE  INPUT 

ALL  VARIABLES  ARE.  FLOATING  POINT  (REAL)  NUMBERS  UNLESS  OTHER¬ 
WISE  SPECIFIED* FLOATING  POINT  NUMBERS  ARE  NUMBERS  WITH  A  DECIMAL 
POINT*  MOST  OF  THE  VARIABLES  WHICH  ARE  FLOATING  POINT  NUMBERS  HAVE 

been  allotted  is  or  is  columns,  if  a  floating  point  number  must 

BE  WRITTEN  IN  ITS  EXPONENTIAL  FORM  (E*G*  ± *XXXXXXXXE±XX  OR  EQUI¬ 
VALENTLY  t  *XXXXXXXX±>XX)  tTHEN  THE  EXPONENT  MUST  BE  RIGHT  MOST  AD¬ 
JUSTED  IN  THE  COLUMNS  ALLOTTED. 

FIXED  POINT  VARIABLES  MUST  BE  RIGHT  MOST  ADJUSTED  IN  THE  COL¬ 
UMNS  SPECIFIED. IF  ON  CARD  2*  NOSAT  =1*  THEN  1  MUST  BE  PUNCHED  IN 
COLUMN  5.  IF  NOSAT  =12 * THEN  1*2  MUST  BE  PUNCHED  IN  COLUMNS  4*5 
RESPECTIVELY. 


B.  INPUT  SYMBOLS 


Caro  variable  definition  columns 

1  TITLE  ANY  DESCRIPTIVE  TITLE  1-66 

2  A)  NOSAT  SATELLITE  NUMBER. MUST  AGREE  WITH  CORRESPONDING  1-5 

COLUMNS  OF  OBERVATION  CARDS.  (FIXED  POINT) 

B)  NMS  NUMBER  OF  STATIONS  INPUT.  liNMS560  (FIXED  POINT)  7-8 


C)  NOB  =1 'REJECT  ANY  OBSERVATION  CARD  WITH  ELEVATION  11 

ANGLE  LESS  THAN  ELEMIN  OR  GREATER  THAN  ELEMAX. 

=0»D0  NOT  REJECT  ANY  OBSERVATION  CARD 

D)  NAT  =1 'CORRECT  ELEVATION  AND/OR  RANGE  AND/OR  RANGE  14 

RATE  FOR  REFRACTION. 

=0»N0  CORRECTION 

E)  MERASE  =1» INPUT  ZONAL  HARMONICS (2-6) »REARTH»F»OMU.  17 

=0»DO  NOT  INPUT  ZONAL  HARM0NlCS(?-6) »REARTH»F#OMU. 

****  SEE  CARDS  11  AND  12  **** 

F)  KOrtSPR  =1 'PRINT  OUT  FILTERING  OUTPUT  20 

=0.00  NOT  PRINT  OUT  FILTERING  OUTPUT 

G)  KCOUNT  NUMBER  OF  OBSERVATION  CARDS  TO  BE  FILTERED  21-23 


MORE  THAN  ONCE.  0 tKCOUNT*r200  IF  KCOUNT  IS 
ZERO  THE  PROGRAM  WILL  FILTER  ALL  OBSERVATION 
CARDS  ONLY  ONCE.  IF  KCOUNT  IS  NOT  EQUAL  TO 
ZERO  THE  FIRST  KCOUNT  CARDS  WILL  BE  FILTERED 
THRICE  ANO  ANY  REMAINING  CARDS  WILL  BE  FILTERED 
ONCE. 

H)  ISMOOH  =1.  APPLY  A  SMOOTHING  TECHNIQUE.  SMOOTHING  WILL  26 

OCCUR  ONLY  IF  THE  EPHEMERIS  TIME  EXTENDS  BEYOND 
THE  TIME  OF  THE  LAST  OBSERVATION  CARD 
=0»  NO  SMOOTHING 

I)  NSTPRT  NUMBER  OF  STATIONS  TO  BE  CONSIDERED  DURING  28-2Q 

EPHEMERIS  PRINTOUT.  OfNSTPRT^NMS 

IF  NSTPRT  EQUALS  ZERO  NO  OBSERVATION  DATA  WILL  BE 

COMPUTED. 

IF  NSTPRT  IS  NOT  ZERO'  THE  FIRST  NSTPRT  STA¬ 
TIONS  INPUT (SEE  CARD  3)  WILL  BE  CONSIDERED  FOR 
OBSERVATION  DATA  COMPUTATION.  NOTE'  IF  KEYTAP 
EQUALS  0  OR  -1.  AND  JTYPRT  EQUAL  3*4  OR  5»  THEN 
NSTPRT  MUST  NOT  EQUAL  ZERO.  (  KEYTAP  AND  JTYPRT 
ARE  DEFINED  ON  CARD  13) 


caro  variable  definition  columns 

J)  PASS  =1.0 »  SINGLE  CORRECTION  OF  RANGE  FOR  REFRACTION.  30-32 

=2*0 '  DOUBLE  CORRECTION  OF  RANGE  FOR  REFRACTION. 

K)  TNEXOS  IF  POSITIVE#  NIGHT  EXOSPHERIC  TEMPERATURE  (DEG  K>  45-59 

IF  NEGAIIVE#  10.7  CM  SOLAR  FLUX.  PROGRAM  WILL 
SET  IT  POSITIVE  AND  CONVERT  IT  TO  NIGHT  EXOSPHE¬ 


RIC  TEMPERATURE. 

L)  IREVU  REVOLUTION  NUMBER  AT  THt  TIME  OF  THE  INPUT  60-65 

ORBITAL  ELEMENTS  (FIXED  POINT) 

M)  TIMTOL  SYSTEM  TIMING  ERROR  (SEC)  66-74 

3  STATION  DATA 

A)  NUMSTA ( I )  STATION  NUMBER.  MUST  AGREE  WITH  CORRESPONDING  2-5 

COLUMNS  OF  OBSERVATION  CARDS  (FIXED  POINT) 

B)  NSo(l)  =1 'ELEVATION  ANO  AZIMUTH  ARE  MEASURED  WITH  8 

RESPECT  TO  GEOCENTRIC  SYSTEM 

=2 'ELEVATION  ANU  AZIMUTH  ARE  MEASURED  WITH 

RESPECT  TO  GEODETIC  SYSTEM 

C)  SLON(I)  STATION  LONGlTUuE  (DEGREES)  POSITIVE  TO  WEST  9-23 

FROM  GREENWICH 

D)  PHILAT(I)  STATION  GEODETIC  LATITUDE  (DEGREES)  24-38 

E)  ALT ( I )  STATION  ALTITUDE  (METERS)  39-53 


***  I=i'...»NMS.  ALL  STATIONS  DEFINED  ON  THE  OBSER¬ 
VATION  CARDS  MUST  bE  REPRESENTED  IN  THIS  SET. 

H  SATELLITE  DATA 

A)  UMSAT  MASS  (KG)  1-15 

B)  SSAT  REFERENCE  AREA  (METERS  SQUARE)  16-30 

C)  cdkag  drag  coefficient  31-45 
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CARD 

VARIABLE 

DEFINITION 

COLUMNS 

h 

DATE  AND 

UNIVERSAL  TIME  OF  THE  INPUT  ORBITAL 

elements 

A) 

MOHS ( 2 ) 

YEAR#LAST  2  DIGITS  OF  19XX  (FIXED 

POINT) 

2-3 

B) 

MOHS  <  3 ) 

MONTH  (FIXtD  POINT) 

5-6 

C) 

MOHS (4) 

DAY  (FIXED  POINT) 

8-9 

0) 

^DAT(1 ) 

HOUR  OF  THE  FORM  XX. U 

11-14 

E) 

^DAT(2> 

MINUTE  OF  THE  FORM  XX. 0 

16-19 

F) 

ZDAT (3) 

SECOND  OF  THE  FORM  XX. XXX 

21-26 

t> 

ORhITAl  ELtMENTS 

A) 

AXSFMI 

SEMI-MAJOR  AXIS  (KILOMETERS) 

1-13 

B) 

tCLEN 

ECCENTRICITY 

14-26 

C) 

oInCl 

INCLINATION  (DtGREES) 

27-39 

D) 

trfASC 

RIGHT  ASCENCION  OF  ASCENDING  NODF 

(DEGREES) 

40-52 

E) 

K»PERi 

ARGUMENT  OF  PERIGEE  (DEGREFS) 

53-65 

F) 

XML  AN 

MEAN  ANOMALY  (DtGREES) 

66-78 

/  iNlT AAl  FSIIMA1E  OF  POSITION  AnO  VELOCITY  ERRORS 


A) 

PMAT(l.l) 

estimated 

ERROR 

IN 

X(0) 

COORDINATE 

(KM) 

1-13 

B) 

PMAT (2*2) 

FSTlMATtD 

ERROR 

IN 

Y  ( 0 ) 

COORDINATE 

(KM) 

14-26 

C) 

PM AT (3»3) 

ESTIMATtD 

error 

IN 

*(U> 

COORDINATE 

(KM) 

27-39 

0) 

PMAT ( 4 » 4 ) 

ESTIMATED 

ERROR 

IN 

X(0> 

COORDINATE 

(KM/SEC) 

40-52 

E) 

PMAT (5*5) 

ESTIMATtD 

EKROR 

IN 

Y  ( 0 ) 

COORDINATE 

( KM/SEC ; 

53-65 

F) 

PMAT (6»6 ) 

ESTIMATtD 

ERROR 

IN 

2(0) 

COORDINATE 

(KM/SEC) 

66-7B 
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CARO  VARIABLE  DEFINITION 


COLUMNS 


8  INITIAL  ESTIMATE  OF  MEASUREMENT  ERRORS  (STANDARD  DEVIATIONS) 

A)  OP(l)  ERROR  IN  DECLINATION  (SECONOS  OF  ARC)  1-8 

B)  QP(  2)  ERROR  IN  RIGHT  ASCENSION  (SECONDS  OF  ARC)  9-16 

C)  0(1)  ERROR  IN  ELEVATION  (SECONDS  OF  ARC)  17-24 

D)  0(2)  ERROR  IN  AZIMUTH  (SECONDS  OF  ARC)  25-32 

E)  0(3)  ERROR  IN  RANGE  (KM)  33-40 

F)  0(4)  ERROR  IN  RANGE  RATE  ( KM/SEC )  41-48 

G)  0(5)  ERROR  IN  ELEVATION  RATE  (SECONDS  OF  ARC/SEC)  49-56 

H)  0(6)  ERROR  IN  AZIMUTH  RATE  (SECONDS  OF  ARC/SEC)  57-64 

I)  O'?)  ERROR  IN  RANGE  ACCELERATION  (KM/SEC*)  65-72 

9  CT ( 1-9)  NUMBER  OF  SIGMAS  CONSIDERED  FOR  REJECTION  1-8 

OF  OBSERVATIONS  DURING  THE  FIRST  FILTERING.  9-16 

CT ( 1-9)  REFER  TO  DECLINATION* RIGHT  ASCENSION*  17-24 

ELEVATI0N*AZIMUTH*RANGE*RANGE  RATE*ELEVATION  25-32 

RATE* AZIMUTH  RATE  ANO  RANGE  ACCELERATION  *  RES-  33-40 

PECTIVELY  41-48 

49-56 

57-64 

65-72 

10  A)  ELEMIN  MINIMUM  ELEVATION  ANGLE  ALLOWED  (DEGREES)  1-15 

B)  ELEMAX  MAXIMUM  ELEVATION  ANGLE  ALLOWED  (DEGREES)  16-30 

***  INPUT  CARO  10  IF  AND  ONLY  IF  NOB  IS  EQUAL  TO  1 

11  ZONHAR (2-6)  COEFFICIENTS  OF  ZONAL  HARMONICS  J(2)  1-15 

J(3)  16-30 

J(4)  31-45 

J( 5)  46-60 

J(6)  61-75 
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CAR 


CARD 

VARIABLE 

DEFINITION 

COLUMNS 

12  A) 

REARTH 

MEAN  EQUATORIAL  EARTH  RADIUS 

(KM) 

1-15 

B) 

F 

FLATTENING  OF  THE  EARTH 

16-30 

C) 

OMU 

EARTH  S  GRAVITATIONAL  NUMBER 

(KM**3/SEC**2) 

3131-45 

*♦**  INPUT 

CARDS  11*12  IF  AND  ONLY  IF  MERASE  IS  EQUAL  TO  1 

IF  CARDS  Hr  12  ARE  NOT  INPUT  THEN  THE  PROGRAM  MILL  SET 
Z0NHAR(2-6)  =• 10827E-2*  -2.4E-6*  -1.6E-6*  -.02E-6*  .7E-6 


RESPECTIVELY*  ANO 

HEARTH  =6378.165*  F  =.00335233*  OMU  =398603.20 

CARDS  13  AND  14  REFER  TO  EPHEMER I S  PRINTOUT 

13  A)  DPRlNT  PRINT  INTERVAL  (SECONDS)  1-15 

*♦*  THE  FOLLOWING  6  VARIABLES  REPRESENT  THE  INITIAL  PRINT  TIME 

(DATE. TIME  OF  OAY)  OF  THE  EPHEMER IS. 

B)  JMOPRT(l)  MONTH  (FIXED  POINT)  17-18 

C)  JDYPRT ( 1 )  DAY  (FIXED  POINT)  20-21 

D)  JYRPRT(l)  YEAR .LAST  2  DIGITS  OF  19XX  (FIXED  POINT)  23-24 

E)  HRPRT(l)  HOUR  OF  THE  FORM  XX.O  26-29 

F)  XMIPHT(l)  MINUTE  OF  THE  FORM  XX.O  31-34 

G)  SECPRT(l)  SECOND  OF  THE  FORM  XX.XXX  36-41 

***  THE  FOLLOWING  6  VARIABLES  REPRESENT  THE  FINAL  PRINT  TIME 

(OATE.T1ME  OF  DAY)  OF  THE  EPHEMER IS. 

H)  JMOPRT ( 2 )  MONTH  (FIXED  POINT)  43-44 

I)  JDYPRT(2)  DAY  (FIXED  POINT)  46-47 

J)  JYRPRT(2)  YEAR .LAST  2  DIGITS  OF  19XX  (FIXED  POINT)  49-50 

K)  HRPRT(2)  HOUR  OF  THE  FORM  XX.O  52-55 
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*R0  VARIABLE  DEFINITION 


COLUMNS 


L)  XMIPRT(2)  MINUTE  OF  THE  FORM  XX. 0  57-60 

M)  SECPRT(2)  SECOND  OF  THE  FORM  XX. XXX  62-67 

***  EXAMPLE  IF  THE  VARIABLES  A)  -  M)  ON  CARD  13  WERE 
60.0  6  21  65  2.0  50.0  0.0  6  21  65  23*0  50.0  10.0 

THIS  WOULD  MEAN  TO  PRINT  EVERY  60  SECONDS  (DPRINT) 

FROM  JUNE  21*1965  AT  2<HR)  50 (MIN)  0(SEC.) 

TO  JUNE  21*1965  AT  23<HR)  50(MIN)  10<SEC.) 

N)  NOSPRI  NUMBER  OF  SPECIAL  PRINT  TIMES  REQUESTED.  69-70 

MUST  BE  LESS  THAN  21  (FIXED  POINT) 

O)  jtyprt  indicates  type  of  ephemeris  print-out  requested  73 

=0*NO  EPHEMERIS.  IN  THIS  CASE  ALL  VARIABLES  ON 
CARD  13  MAY  BE  INPUT  AS  ZERO. 

=1*PRINT  POSITION  AND  VELOCITY 
=2* PRINT  POS I T I ON  *  VELOC I T Y  AND  OSCULATING 
ELEMENTS 

=3* PRINT  POSITION*VELOCITY*OSCULATING  ELEMENTS 
AND  STATION  OBSERVATION  OATA. 

=4* PRINT  POSIT ION* VELOCITY  AND  STATION  OBSERVA¬ 
TION  DATA. 

=5*PRINT  STATION  OBSERVATION  DATA 

P)  KEYTAP  INDICATES  TYPE  OF  EPHEMERIS  TAPE  REQUIRED  75-76 

=-l*BCD  TAPE  ONLY 
=  0* BINARY  AND  BCD  TAPE 
=  1*BINARY  TAPE  ONLY 

THE  BCD  TAPE  IS  WRITTEN  ON  U02* 

THE  BINARY  TAPE  IS  WRITTEN  ON  U06(USE  800  BPI  DENSITY  TAPE) 

IF  KEYTAP  =0  OR  1*  THEN  THE  BINARY  TAPE  WILL  CONTAIN  ALL  THE 
POSSIBLE  DATA  AT  ANY  PRINT  TIME  I.E.  IT  ASSUMES  JTYPRT  =3 


Caro  variable  definition  columns 

14  SPECIAL  PRINT  TIMES  (DATE* TIME  OF  DAY)  INPUT  IF  NOSPRI  IS 
GREATER  THAN  ZERO. 

A)  JSPMO(i)  MONTH  (FIXED  POINT)  2-3 

B)  JSPDA(l)  DAY  (FIXED  POINT)  5-6 

C)  JSPYK(I)  YEAR 'LAST  2  DIGITS  OF  19XX  (FIXED  POINT)  6-9 

D)  SPHR(I)  HOUR  OF  THE  FORM  XX .0  11-14 

E)  SPMI(I)  MINUTE  OF  THE  FORM  XX .0  16-19 

F)  SPSEC(l)  SECOND  OF  THE  FORM  XX.XXX  21-26 

«**  1=1 » ... 'NOSPRI 

THESE  SETS  OF  SPECIAL  PRINT  TIME  MUST  FALL  WITHIN  THE 
TIME  INTERVAL  OEFINEO  BY  VARIABLES  B)  -  M)  ON  CARD  13 

***  EXAMPLE'  IF  UN  THE  ABOVE  EXAMPLE  NOSPRU1  AND  CARD  14 


CONTAINS  6  21  65  5.0  6.0  12.398'THEN  A  PRINTOUT  WOULD  ALSO 

OCCUR  FOR  JUNE  21'1965  AT  5(HR)  6(MIN)  12.398(SEC) 

lb  OBSERVATION  CARDS. MUST  BE  INPUT  IN  INCREASING  ORDER  WITH 
RESPECT  TO  TIME.  CARDS  OUT  OF  SEQUENCE  WILL  BE  REJECTED. 
DECIMAL  POINTS  ARE  NOT  PUNCHED  ON  THESE  CARDS. THEY  ARE 
IMPLIED  BY  THE  FORTRAN  FORMAT. 


A)  MERASE(l) 

SATELLITE  NUMBER 

2-6 

B)  MOBS ( 1 ) 

STATION  NUMBER 

7-9 

C)  MOBS ( 2 ) 

LAST  2  DIGITS  OF  YEAR  19XX 

10-11 

D)  DAYS 

DAY  OF  YEAR 

12-14 

E)  ZOAT(l) 

HOUR  OF  DAY 

15-16 

F)  ZD AT (2) 

minute  OF  HOUR 

17-16 

Caro  variable  definition  columns 

G)  *DAT(3)  SECOND  OF  MINUTE.  DECIMAL  POINT  IMPLIED  19-23 

BETWEEN  COLUMNS  20  AND  21 

H)  MERASE15)  FIRST  DIGIT  OF  ELEVATION  OR  DECLINATION. (DEG)  24 

USE  AN  OVERPUNCH  (11  PUNCH)  FOR  NEGATIVE 
QUANTITIES 

I)  OBSM(l)  REMAINING  DIGITS  OF  ELEVATION  OR  DECLINATION.  25-29 

DECIMAL  POINT  IMPLIED  BETWEEN  COLUMN  25  AND 
26. 

J)  ERASE ( I )  FIRST  2  DIGITS  OF  AZIMUTH  (DEG) »  OR  31-32 

HOUR  OF  RIGHT  ASCENSION 

K)  ERASE (2)  NEXT  2  DIGITS  OF  AZIMUTH. DECIMAL  POINT  IMPLIED  33-34 

BETWEEN  COLUMNS  33  AND  34*  OR 
MINUTE  OF  RIGHT  ASCENSION 

L)  ERASE(3)  LAST  3  DIGITS  OF  ELEVATION'OR  35-37 

SECONDS  OF  RIGHT  ASCENSION. DECIMAL  POINT 
IMPLIED  BETWEEN  COLUMNS  36  AND  37 

M)  OBSM13)  RANGE  (KM)  39-45 

N)  1EX  EXPONENT  46 


AN  EXPONENT  OF  ZERO  (0)  IMPLIES  THE  DECIMAL 
POINT  BETWEEN  COLUMNS  40  AND  41 'AN  EXPONENT 
OF  ONE  (1)  IMPLIES  THE  DECIMAL  POINT  BETWEEN 
COLUMNS  41  AND  42  ETC.  THE  MAXIMUM  ALLOWABLE 
EXPONENT  IS  FIVE  (5) 


O)  0BSM(4)  FIRST  DIGIT  OF  RANGE  RATE  (KM/SEC).  USE  AN  48 

OVEkPUNCH  (11  PUNCH)  FOR  NEGATIVE  QUANTITIES 

P)  ERASE (b)  REMAINING  DIGITS  OF  RANGE  RATE. DECIMAL  POINT  49-54 

IMPLIED  BETWEEN  COLUMNS  49  AND  50 

Q)  QBSM(5>  FIRST  DIGIT  OF  ELEVATION  RATE  (DEG/SEC).  USE  56 

AN  OVERPUNCh  FOR  negative  quantities 

R)  ERASE (e>)  REMAINING  DIGITS  OF  ELEVATION  RATE. DECIMAL  57-60 

POINT  IMPLIED  BETWEEN  COLUMNS  5b  AND  57 
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CARO  VARIABLE  DEFINITION 


COLUMNS 


S)  OBSM ( 6 )  FIRST  DIGIT  OF  AZIMUTH  RATE  (OEG/SEC)  USE  62 

AN  OVERPUNCH  FOR  NEGATIVE  QUANTITIES 

T)  ERASE (7)  REMAINING  DIGITS  OF  AZIMUTH  RATE. DECIMAL  63-66 

POINT  IMPLIED  BETWEEN  COLUMNS  62  AND  63 

U)  OBSM ( 7 )  FIRST  DIGIT  OF  RANGE  ACCELERATION  (KM/SEC2)  68 

USE  AN  OVERPUNCH  FOR  NEGATIVE  QUANTITIES 

V)  ERASE (8)  REMAINING  OIGITS  OF  RANGE  ACCELERATION • DEC I -  69-72 

MAL  POINT  IMPLIED  COLUMNS  67  AND  68 

W)  OBSNO  CODE  WHICH  DESIGNATES  TYPE  OF  SIMULTANEOUS  75-76 

OBSERVATIONS  IS  IN  COLUMN  75. THIS  CODE  WILL 
BE  REDEFINED  ANO  STORED  IN  THE  VARIABLE  NTP. 

(SEE  NTP  IN  LIST  OF  SYMBOLS) 

COL. 75  TYPES  OF  OBSERVATIONS  NTP 

=1  ELEVATION  AND  AZIMUTH  2 

=2  ELEVATION* AZIMUTH  AND  RANGE  3 

=3  ELEVATION* AZIMUTH* RANGE  AND  4 

RANGE  RATE 

—4  ELEVATION* AZIMUTH* RANGE* RANGE  RATE*  7 

ELEVATION  RATE* AZIMUTH  RATE  AND 
RANGE  ACCELERATION 

=5  DECLINATION  AND  RIGHT  ASCENSION  1 


THE  PROGRAM  CONVERTS  THE  RIGHT  ASCENSION 
TO  RADIANS  BY  THE  EQUATIONS 

OBSM(2)= (ERASE (1)  +  ERASE(2)/60.0  ♦ 

ERASE(3)/3600.0) *15.0*. 017453292 


COOL  OF  THE  CELESTIAL  REFERENCE  SYSTEM (YEAR 
OF  EQUINOX)  IS  IN  COLUMN  76*  THIS  IS  USED  ONLY 
WHEN  THt  DECLINATION  OR  RIGHT  ASCENSION  ARE 
INPUT .THE  COOE  IS  AS  FOLLOWS 

COL .76  YEAR 

=0  YEAR  OF  DATE  (M0BS(2>) 

=1  1900.0 

=2  1920.0 

=3  1950.0 

=4  1975.0 

=5  2000.0 

=b  1850.0 

=7  1855.0 

=8  1875.0 

=9  1960.0 

this  code  number  will  be  stored  in  the  variable 

NtQ  (SEE  NEQ  IN  LIST  OF  SYMBOLS) 

lb  A  CARD  WITH  THE  WORD  ENDDAT  IN  COLUMNS  73-78.  THIS  CARD  73-78 

FOLLOWS  THt  LAST  OBSERVATION  CARD 

17  A  CARD  WITH  THt  PHRASE  SEND  (IN  COLUMNS  1-3)  OF  (IN  COLUMNS  1-14 

b-b)  PNOBLtM  (IN  COLUMNS  8-14) 


20 


IV.  PROGRAM  OUTPUT 


A. OUTPUT  OF  THt  INPUTCEXCEPT  FOR  THE  OBSERVATION  CARDS) 

THE  OUTPUT  FOR  THIS  SECTION  IS  DONE  IN  MAIN  PROGRAM  AFILT2* 


INITIAL  CONDITIONS 

REVOLUTION  NUMBER-IREVO 

DATE  AND  TIME 

MOBS ( 2 )  MOBS (3)  M0BS(4)  ZDAT(l)  ZDAT(2)  ZDAT(3) 

ORBITAL  ELEMENTS 

AXSEMI  ECCEN  OlNCL  MASC  tfPERI  XMEAN 

POSITION  AND  VELOCITY  ERRORS 

PMAT ( 1 » 1 )  PMAT (2»2)  PMAT(3>3)  PMAT(4»4)  PMAT(5#5)  PMAT(6»6) 

MEASUREMENT  ERRORS 

QP(1)  QP ( 2)  0(1)  Q (2)  0(3)  Q (4) 

0(5)  Q(b)  0(7) 

NUMBER  OF  SIGMAS  TOLERATED  CT(l-9) 

SATELLITE  DATA  EARTH  DATA 

OMSAT  SSAT  CORAG  HEARTH  F  OMU 

ZONAl.  HARMONICS (2—6)  Z0NHAR(l-5) 

PROGRAM  CONTROLS  NOB  NAT  MERASE ( 1 )  KOBSPR  KCOUNT  ISMOOH  NSTPRT 

PASS  TIMTOL 

STATION  DATA 

NUMSTA  NSG  SLON  SLAT  PHILAT  ALT  SRAD 

EPHEMERIS  PRINTOUT  DATA 

JMOPRT(I)  JOAPRT(l)  JYRPRT(l)  HRPRT(l)  XMIPRT(l)  SECPRT(l) 

JM0PRT(2)  JUAPR1 (2)  UYRPRT(2)  HRPRT ( 2 )  XMIPRT(2)  SECPRT(2) 

PRINTOUT  CODE-  JTYPRT  TAPE  INDICATOR-  KEYTAP  PRINT  INTERVAL-  DPRlNT 

SPECIAL  PRINTOUT  TIMES  (IF  ANY) 

JSPMU  JSPDA  JSPYR  SPHR  SPMI  SPSEC 

INITIAL  COMPUTATIONS 

PERIOD  ERASE ( 1 )  TNEXOS 

PNODaL 
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REMARKS  1)  THE  DIMENSION  OF  ALT *STATI0N  ALTITUDE* IS  IN  KILOMETERS ( IT 

is  input  in  meters) 


2)  THE  variables  slat*srao*period  and  pnodal  are  computed. 

(SEE  THE  LISTING  OF  SYMBOLS  FOR  THEIR  DEFINITIONS) 

3)  ERASE <t)  IS  THE  COMPUTED  ECCENTRICITY  (DEFINED  IN  SUBROU¬ 
TINE  KEPLER) 

4)  THE  REMAINING  VARIABLES  ARE  DEFINED  IN  THE  INPUT  LISTING. 


5)  if  any  errors  in  the  input* the  program  will  printout  an 

APPROPRIATE  MESSAGE  AND  CONTINUE  TO  THE  NEXT  PROBLEM. 


O.F iLTtRING  OUTPUT 

THE  FILTERING  OUTPUT  IS  DONE  IN  MAIN  PROGRAM  BFILT2.EACH  PAGE  IS  PRO¬ 
PERLY  LABELED  AND  THE  PHRASE  FORWARD  FILTERING  (OR  BACKWARD  FILTER¬ 
ING)  IS  WRITTEN  ON  IT. 


DJULN  MOBS (2)  MOBS (3)  M0HSC4)  ZUAT(l)  /DAT ( 2 ) 


DV  ( 1 ) 
ERASt ( l ) 
ERASE ( 7 ) 


Dv  (2) 
ERASE (2) 
ERASE («) 


ELEVmTAON  OR 
DECLINATION 
AZIMUTH  OR 


DV  (3) 
ERASE (3) 
ERASE (H) 
MEASURED 
OHSM(l) 


DVP(l) 
ERASE (4 ) 
ERASE ( 10) 


DVP ( 2 ) 
FRA5E (5) 
FRASE(il) 


COMPUTED 
OBSC ( 1 ) 


uRSM ( 2 )  OBSC (2) 


DIFFERENCE 

DIFR(l) 

DIFR ( 2 ) 


7DAT  (3)  NIJMSTA 
DVP<3) 

ERASE (6) 
ERASE( 12) 
TOLERANCE 
TOLER ( 1 ) 

TOLER (2) 


RlGHl 

ascension 

RANGt 

UBSM 

RANGt 

RATE 

URSM 

elevation  RATE 

URSM 

AZIMUTH  RATE 

URSM 

range 

ACCELERATION 

URSM 

(3) 

OBSC (3) 

DIFR (3) 

(4) 

OBSC (4) 

DIFH ( 4 ) 

(5) 

OBSC (5) 

DIFR (5) 

(6) 

OBSC (6) 

DIFR(6) 

(7) 

OBSC (7) 

DIFR  ( 7 ) 

TOLER (3) 
TOLER (4) 
TOLER (S) 
TOLER (6) 
TOLER (7) 


REMAkKS  1)  DJULN  IS  THE  MODIFIED  JULIAN  DATE  OF  THE  OBSERVATION. 

MOBS (2-4)  REPRESENT  THE  DATE  OF  THE  OBSERVATION. 

7DaT ( 1-3)  REPRESENT  THE  TIME  OF  DAY  OF  THE  OBSERVATION. 
NUMSTA  IS  THE  STmTION  NUMBER. 


22 


2)  DV ( 1-3)  ARt  THt  X*Y*Z  COORDINATES  OF  THE  UPDATED  POSITION 
VECToR.  UN  KILOMETERS) 

UVPU-3)  ARE  THE  XfYii'  COORDINATES  OF  THE  UPDATED  VELOCITY 
VECTOR. (1N  KlLOMtTERS/SECOND) 

3)  tRASt(l-o)  ARE  THE  CORRECTIONS  TO  THt  POSITION  AND  VELO¬ 
CITY  VtCl OKS. ( IN  KILOMETERS  AND  KlLOMETERS/SECONO»RES- 
HECTaVELT) 

ERASt ( 7—12)  ARE  THE  STANDARD  DEVIATIONS  OF  THE  POSITION 
ANU  VELOCITY  EKRuRS. ( IN  KILOMETERS  AND  KILOMETERS/SECOND* 
RESPECTI vELY ) 

eRaSc ( 1-12)  ARE  DEFINED  IN  BFILT2* 

4)  THe  mNGUlAr  DATA  OF  THE  OBSERVATIONS  ARE  OUTPUT  IN  DEG¬ 
REES  OR  dEoReES/SECOND. 

THE  RANGE  IS  In  KILOMETERS. 

THt  kANGE  RATE  IS  IN  KILOMETERS/SECOND. 

THE  kANGc  ACCELERATION  IS  IN  KILOMETERS/SECOND2. 

5)  IF  Aim  OBSERVATION  CARD  IS  REJECTED*  THE  FOLLOWING  MESSAGE 
wltL  BE  PRINTED  uUT: 

FOLLOW! Nb  LARD  DiSCAhDtD.ObStRVATlON  DEVIATES  FROM  THEORETICAL 
VALUE  dEYOND  AtLOWABcE  LIMIT • 

THtN  THE  PROGRAM  PRINTS  THE  SAME  INFORMATION  AS  ABOVE » 

except  fur  the  variables  erase ( 1-12 > 

6)  iF  TmE  PkOGRaM  I1EHATES*THEN  THE  PROGRAM  WILL  PRINTOUT: 


ITERATION  CHtCK  ON  DtLTA  POSITION  AND  VELOCITY  FAILED. 

P0S=UV<1)  DV  U)  DV  ( 3)  VEL=DVPU)  UVH2)  DVP(3) 

DEl=eRASe(1)  ERASE(2)  ERASE ( o )  DEL=ErASE(4)  ERASE (5)  EHASE(6) 

posit  ion* velocity  *  and  cokrlction  for  previous  time 

UVU)  DV  (2)  uV  ( 3)  DVP(l)  DVP(?>  DvP(3) 

RASE ( 1 )  RAbE (2)  «AbE(3)  RASE<4)  PASt(5)  RASE(o) 

RAbE l 1-6)  ARE  THE  CORRECTIONS  TO  THE  PREVIOUS  POSITION 
AND  VELOCITY  VECTOR • ( IN  KILOMETERS  AND  KlLOMETERS/SECOND) 

7)  AT  THE  END  OF  tACH  FILTERING  PHASE  THE  PROGRAM  PRINTS: 

covariance  matrix 

PM  AT 1 1  *  J)  1-1  *  ...  *6  Jil  *  .  •  .  *6 
NEW  MEASuREMtNT  ERRORS 
ERASt (1-9) 

NUMBER  OF  SIoMaS  TOLERATED  LT ( 1-9) 
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PMAT  IS  THE  COVARIANCE  MATRIX. 

ERASE ( 1-9)  ARE  THE  IMPROVED  STANDARD  OBSERVATION  ERRORS. 
(SEE  EQUATION  A. 15) 

CT ( 1-9)  ARE  THE  RECOMPUTED  NUMBER  OF  SIGMAS. (SEE  EQUATION 
A.16M) 

THE  PROGRAM  ALSO  PRINTS  THE  COMPUTING  TIME  FOR  THIS  PHASE 
IN  MINUTES 

TIME  FOR  THIS  PHASE=  ERASE ( 1 )  MIN. 

8)  IF  THERE  IS  OVERFLOW  OR  A  DIVISION  BY  ZERO* THE  PROGRAM 
PRINTS  AN  APPROPRIATE  MESSAGE  * DUMPS  CORE* AND  GOES  TO  THE 
NEXT  PROBLEM. 

9)  THE  DIFFERENTIAL  EQUATION  SUBROUTINE  PRINTS  THE  FOLLOW¬ 
ING  ERROR  MESSAGE 

DENS IT  SUBROUTINE  ERROR  NO.RHO.  PROGRAM  WILL  GO  TO  NEXT  PROBLEM. 

IF  RHO  =-1.0* THE  EXOSPHERIC  TEMPERATURE  IS  OUTSIDE 
OF  THE  500° K  -  2400° K  RANGE. 

IF  RHO  =  0.0* THE  ALTITUDE  IS  LESS  THAN  120.0 (KM). 

10)  IF  THE  DIAGONAL  ELEMENTS  OF  THE  COVARIANCE  MATRIX  ARE 
NEGATIVE* THEN  THE  PROGRAM  PRINTS 
COVARIANCE  MATRIX  BECAME  NEGATIVE. 

IF  THIS  NUMERICAL  PROBLEM  REOCCURS  AFTER  DI AGONALIZATION* 
THEN  THE  PROGRAM  WILL  DUMP  CORE* AND  GO  TO  THE  NEXT 
PROBLEM. 

(SEE  EQUATION  51  * REFERENCE  1  ANO  SUBROUTINE  DIACHK) 


11)  IF  THE  SMOOTHING  OPTION  HAS  BEEN  REQUESTED* THE  PROGRAM 
WILL  PRINT  ON  THE  SYSTEM  OUTPUT  TAPE 
SMOOTHING  DATA.  POINT  TO  BE  SMOOTHED  AT  OLDDAY*OLDTIM 

EACH  TIME  THE  PROGRAM  COMPUTES  THE  SMOOTHING  EQUATIONS 
THE  PROGRAM  PRINTS 

MOD. JUL. DATE  =  DJULN  TIME  =  TNORMN(SEC) 

DIFFERENCES  BETWEEN  STORED  AND  INTEGRATED  VALUES  IN  POSI1 ION  AND  VELOCITY 
ERASE ( 1 )  ERASE (2)  ERASE (3)  ERASE (4)  ERASE(5)  ERASE(6) 

CORRECTIONS  IN  POSITION  AND  VELOCITY  AT  THE  INITIAL  DATE  ( I .E.OLDDAY»OLDTIM) 
ERASE ( 1 )  ERASE (2)  ERASE (3)  ERASE (4)  ERASE (5)  ERASE(G) 

NEW  POSITION  AND  VELOCITY  AT  THE  INITIAL  DATE 

DV ( 1 )  DV(2)  DV ( 3)  DVP(l)  DVP(2)  DVP(3) 

OLDDAY  IS  THE  MODIFIED  JULIAN  DATE  OF  THE  LAST  OBSERVATION 

CARD  WHICH  HAS  BEEN  STORED  FOR  SMOOTHING 

OLDTIM  IS  THE  TIME  OF  DAY  IN  SECONDS  CORRESPONDING  TO 

OLDDAY. 

THIS  IS  PRINTED  IN  SUBROUTINE  ORBINT 
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C.EPHEMERIS  OUTPUT 


THIS  OUTPUT  IS  WRITTEN  ON  A  BCD  TAPE (WHEN  REQUESTED  BY  CUSTOMER)  AND 
IS  PRINTED  OUT  IN  SUBROUTINE  PROUT.EACH  PAGE  IS  PROPERLY  LABELED  AND 
THE  WORD  EPHEMERIS  IS  WRITTEN  ON  IT. 


DAYJl  KMOOUT  KOAOUT  KYROUT  KHROUT  KM i OUT  SECOUT  TNORMO  IREV 

THE  REMAINDER  OF  THE  OUTPUT  IS  A  FUNCTION  OF  THE  VARIABLE  JTYPRT. (SEE 
1NPU1  LISTING) 

POSITION  AND  VELOCITY  OUTPUT 


DV  ( 1 ) 

DV(2) 

DV(  3) 

DVP(l) 

DVP(2) 

DVP(3) 

ALTIO 

RADPRT  VTCTPR  GEOCEN 

GEODE T  I 

olamo  msthr 

MSTMIN  STSECO 

OSCULATING  ELEMENTS 

OUTPUT 

AXIMaJ 

ECCENT 

OINCLI 

ASCNOD 

PERIGE 

OMEANA 

station 

OBSERVATION 

DATA  OUTPUT 

NUMSTA 

elrate 

AZRATE 

RARATE 

dcrate 

ELEVAT 

AZIMUT 

RANGES 

RANRAT 

RTASC 

DECLIN 

REMARKS  1)  THE  BINARY  TAPE  FORMAT  PRODUCED  BY  THE  PROGRAM  IS  AS 
FOLLOWS 


FIRS!  RECORD 
WORD 

1.  KF7 

2.  NOSAT 

3.  TiMSEUl) 

4.  DoUPRT ( 1 ) 
b.  TaMSEC(2> 
b.  DuUPRT (2) 
7.  DPRINT 

6.  NOSPR1 


IDENTIFICATION  RECORD 
DEFINITION 

7 • NUMBER  OF  WORDS  REMAINING  IN  THIS  RECORD 
SATELLITE  NUMBER 

TIME  OF  DAY  OF  INITIAL  PRINT  TIME  (SEC) 
MODIFIED  JULIAN  DATE  OF  FINAL  PRINT  TIME 
TIME  OF  DAY  OF  FINAL  PRINT  TIME  (SEC) 

modified  Julian  date  of  final  print  time 

PRINT  INTERVAL  (SEC) 

NUMBER  OF  SPECIAL  PRINT  TIMES 
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SECOND  RECORD 


DATA  RECORD 


WORD 

1.  MERASE(l) 

2.  DmYJL 

3.  KMOOUT 

4.  KUAOUT 
b.  KtROUT 

6.  KHROUT 

7.  KMlOUT 
b.  StCOUT 
9.  TNORMO 

10.  DV  ( 1 ) 

11.  DV(2) 

12 .  DV(3) 

13.  DvP(l) 

14.  DvP(2) 
lb.  DVP(3) 
lb.  AlTIO 
17.  RaDPRT 

lti.  VTOTPK 

19.  GtOCEN 

20.  GtODET 

21.  OlAMO 

22.  MbTHR 

23.  MSTMIN 

24.  S1SEC0 
2b.  AxIMAJ 

26.  EcCENT 

27.  OlNCLl 
2b.  ASCNOD 

29.  PERIGE 

30.  ONE ANA 

31.  IKEV 


DEFINITION 

NUMBER  OF'  kkOKDS  REMAINING  IN  THIS  RECORD 
=32  ♦  ll*lJ 

MODIFIED  JULIAN  DATE  OF  THE  EPHEMERIS  PRINT  TIME 
CALENDAR  MONTH 
CALENDAR  DAY 

CALENDAR  YEAR  (LAST  2  DIGITS  OF  19XX) 

HOUR  OF  DAY 
MInUTE  OF  HOUR 
SECONDS  OF  MINUTE 

TIME  OF  DAY  IN  SECONDS  CORRESPONDING  TO  DAYJL 
X  COORDINATE  OF  POSITION  VECTOR  (KM) 

Y  COORDINATE  OF  POSITION  VECTOR  (KM) 

2  COORDINATE  OF  POSITION  VECTOR  (KM) 

X  COORDINATE  OF  VELOCITY  VECTOR  (KM/SEC) 

Y  COORDINATE  OF  VELOCITY  VECTOR  (KM/SEC) 

Z  COORDINATE  OF  VELOCITY  VECTOR  (KM/SEC) 
SATELLITE  ALTITUDE  (KM) 

DISTANCE  OF  THE  SATELLITE  FROM  THE  CENTER 
OF  THE  EARTH  (kM) 

VELOCITY  (KM/SEC) 

GEOCENTRIC  LATITUDE  (DEG) 

geudetic  latitude  (deg) 

SATELLITE  LONGITUDE  (DEG) 

HOUR  OF  GREENWICH  MEAN  SIDEREAL  TIME 
MINUTE  OF  GREENWICH  MEAN  SIDEREAL  TIME 
SECONDS  OF  GREENWICH  MEAN  SIDEREAL  TIME 
SEMI-MAJOR  AXIS  (KM) 

ECCENTRICITY 
INCLINATION  (DEG) 

RIGHT  ASCENSION  OF  ASCENDING  NODE  (DEG) 

ARGUMENT  OF  PERIGE  (DEG) 

MEAN  ANOMALY  (UEG) 
revolution  NUMBER 
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32. 

Io 

NUMBER  OK  STATIONS  IN  THIS  DATA  RECORD  OBS¬ 
ERVING  THE  SATELLITE  (IJilO) 

33. 

HYPE 

=i»final  data  record  for  this  EPHEMERIS  PRINT 

TlMt 

=-l# ANOTHER  DATA  RECORD  TO  FOLLOW  FOR  THIS 
EPHEMERIS  PRINT  TIME 

34. 

NUMSTA 

STATION  NUMBER 

35. 

ELRATt 

ELEVATION  KATE  (DEG/SEC) 

3o . 

AERATE 

AZIMUTH  RATE  (DEG/SEC) 

37. 

RaRATE 

RIGHT  ASCENSION  RATE  (DEG/SEC) 

3b . 

DCRATt 

DECLINATION  RATE  (DEG/SEC) 

39. 

EcEVAT 

ELEVATION  (DEG) 

40. 

A^IMUT 

AZIMUTH  (DtG) 

41. 

ranges 

RANGE  (KM) 

42. 

RaNRAT 

RAimGE  KATE  (KM/SEC) 

43. 

R1ASC 

RIGHT  ASCENSION  (DEG) 

44. 

OeClIn 

DECLINATION  (DEG) 

1HE  SEQUENCE  OF  WORDS  34  TO  44  IS  REPEATtD  IJ  TIMES. 


THIRu  RECORD  DmTA  RECORD  WHICH  CONTAINS  STATION  OBSERVATION 

DATA  ONLY. THIS  RECORD  IS  WRITTEN  WHENEVER  THE  FIRST 
DATA  RECORD  DoES  NOT  CONTAIN  ALL  OF  THE  station 
OdSERvATIUN  OATA. 

WORD  DEFINITION 

1.  MtRASt(l)  NUMBER  OF  WORDS  REMAINING  IN  THIS  RECORD 

=2  ♦  11*1J 

2.  I o  SEE  WORD  32  ABOVE 

3.  I1YPE  SEt  WORD  33  ABOVE 

4.  -  14.  SEt  WORDS  34  -  44  ABOVE 


I  HE  SEQUENCE  OF  WORDS  4  TO  14  IS  RtPEATEU  IJ  TIMES. IF  ITYPE  IS 
STILL  tOUAL  TO  -1 # THtN  THIS  DATA  RECORD  IS  REPEATED. 


1  HE  SECOND  RECORD  (AND  THE  THIRD  RECORD  WHENEVER  NECtSSARY) 
IS  RtPtATEU  FOR  tVERY  PRINT  TIME. 


FOURTH  RECORD 

IDENTIFICATION 

RECORD  FOR  THE  END  OF  THE  PROBLEM 

WORD 

DEFINITION 

1.  KF1 

1 

2.  BLANKS 

OCTAL  NUMBER 

606060606060 

AN  END  OF  FILE  FOLLOWS  THIS  RECORD. 


2)  IF  A  BINARY  TAPE  HAS  BEEN  REQUESTED* THE  PROGRAM  WILL 
PRINT  ON  THE  SYSTEM  OUTPUT  TAPE 

AX1MAJ  tCCEN  OINCLI  ASCNOD  PERIGE  OMEANA 
AREV 

THESE  ARE  THE  OSCULATING  ELEMENTS  AND  REVOLUTION  NUMBER 
FOR  THE  FINAL  PRINT  TIME. 

3)  IF  ANY  ERROR  MESSAGES*  THE  EPHEMERIS  OUTPUT  WILL  NOT  BE 
COMPLETE. 

4)  IF  THEHE  IS  MORE  THAN  ONE  FILTERING  AND  ALL  OBSERVATION 
CARDS  WERE  REJECTED* THEN  THERE  IS  NO  EPHEMERIS  COMPUTA¬ 
TION. 
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*  A  A  o  # 
0*4*  A  O 

•  *  4  •-*  X  A  • 

,N  W  4  O  J  N  m 

tx  4  *  A  4  4  X 

iX  A  fM  »  4  A  X 

|4  •  •  •  4  •  4 

0  0  0  4  0 

tit 


X* 


0« 106695J&C  C«.  0.6400j710fe  04  -0.j7«tt7l66£  04  -0.7»3*»0«3E  01  -0.34147414E  00  -0.999M1954E  00 

0*1 1411 13JE  C4  0» 731 30 J13t  04  0.76217970E  01  -30.2*0449  -  30.422720  222.4097*  20  11  47. *21 

0.830304426  C4  0.164427*26  00  G.3247S039E  02  0.195V2344E  03  0.3110*905E  03  0.31300290E  03 

STATIUM  Nu.  •  33  -0.19*082876  00  0.277189296  01  0.394333376  00  -0.72*40*476-01 

<r.FJ74*Jl96  C 2  0.'47C4a339r  FI  57TT4U95IIE  04 - =0767225 17TE "00 - 0:841472982  02  - -0:263122311  02 


V.  LIST  OF  SYMBOLS 


THE  FOLLOWING  IS  A  LIST  OF  SYMBOLS. THE  METRIC  SYSTEM  OF  UNITS 
IS  USEO  IN  THE  PROGR AM— K ILOGR AMS  *  K I LOMETERS * METERS  AND  SECONDS. 

ALL  ANGULAR  QUANTITIES  ARE  IN  RADIANS  DURING  THE  COMPUTATIONAL  PRO¬ 
CESS  AND  THEY  ARE  INPUT  OR  OUTPUT  IN  DEGREES  OR  SECONDS  OF  ARC. THE  CON¬ 
VENTIONS  USED  IN  THE  DEFINITIONS  ARE 

1)  INTERMEDIATE  VARIABLE— THE  VALUE  OF  AN  EQUATION (OR  PART  OF  AN 
EQUATION)  USED  FOR  INTERMEDIATE  CALCULATION. 

2)  ASSIGNED  VARIABLE— THE  VARIABLE  HAS  SEVERAL  DEFINITIONS. IT  IS 
PROPERLY  DEFINED  IN  THE  MAIN  PROGRAMS  OR  SUBROUTINES  WHERE  IT 
IS  USED. 

3)  SEE  INPUT  LISTING— INPUT  VARIABLE  DEFINED  IN  THE  INPUT  LISTING. 

4)  MODIFIED  JULIAN  DATE— THE  NUMBER  OF  INTEGRAL  DAYS  SINCE  JANUARY  1' 
1950  (OhUT). 

5)  EQUA — EQUATION.  IF  NO  REFERENCE  IS  GIVEN* THE  EQUATION  WILL  BE 
FOUND  IN  THIS  REPORT. 

6)  REF— REFERENCE. 
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VARIABLE 

EGUA 

REF 

DEFINITION 

A  ( 1 ) 

=1.0»HIGHER  ORDER  TERMS  OF  THE  EARTh's  GRA¬ 
VITATIONAL  POTENTIAL  FUNCTION  ARE  INCLUDED 

IN  THE  MATHEMATICAL  MODEL. 

=0.0*  NOT  INCLUDED 

A (2-7)  A.33-A.3to 

INTERMEDIATE  VARIABLE 

A1 

INTERMEDIATE  VARIABLE 

A1M1 

5 

2 

INTERMEDIATE  VARIABLE 

A1P1 

5 

2 

INTERMEDIATE  VARIABLE 

A5 

4 

2 

INTERMEDIATE  VARIABLE 

AELIP 

A. 2b 

INTERMEDIATE  VARIABLE 

AH 

INTERMEDIATE  VARIABLE 

ALP 

INTERMEDIATE  VARIABLE 

ALT ( I ) 

1=1 *  •  .  . »NMS 

5-7 

1 

SEE  INPUT  LISTING.  THE  PROGRAM  CONVERTS 

ALT ( I )  TO  EARTH  RADII. 

ALT  I 

ASSIGNED  VARIABLE 

ALT  10 

A. 58 

SATELLITE  ALTITUDE  (KM)rOUTPUT 

ANGDAT (111) 

STORAGE  ARRAY  FOR  COMPUTED  STATION  OBSERVATION 
DATA 

ARG 

27 

1 

REFRACTION  CORRECTION  FACTOR  FOR  RANGE  RATE 

ASCNOD 

A.  38 

RIGHT  ASCENSION  OF  THE  ASCENDING  NODE (DEG) » 
OUTPUT 

AST I ME 

A. 54 

INTERMEDIATE  VARIABLE 

AX1MAJ 

A. 32 

SEMI-MAJOR  AXIS  (KM)*OUTPUT 

AXSEMI 

SEE  INPUT  LISTING 

AZIMUT 


151  1 


AZIMUTH  (DEG) » OUTPUT 


VARIABLE 

EOUA 

REF 

DEFINITION 

AZKATE 

13b 

1 

AZIMUTH  RATE  (DEG/SEC ) t OUTPUT 

b 

M 

2 

-•78539816  RADIANS 

bl 

b 

2 

INTERMEDIATE  VARIABLE 

63 

b 

2 

INTERMEDIATE  VARIABLE 

6ELIP 

A. 29 

INTERMEDIATE  VARIABLE 

bJK ( 1-3) 

STORAGE  CELLS 

bLANKS 

THE  OCTAL  NUMBER  606060606060 

BLOCK (br200) 

STORAGE  CELLS  FOR  X*Y»Z»X»Y»Z  (FOR  200  CON 
SECUTIVE  PRINT  TIMES) 

bN 

A.23A 

INTERMEDIATE  VARIABLE 

bNO 

A. 22 

APPROXIMATION  TO  THE  NUMBER  OF  REVOLUTIONS 
SINCE  THE  FIRST  EQUATORIAL  CROSSING 

bNl 

A  .23b 

INTERMEDIATE  VARIABLE 

bOX(brS) 

A.16A 

STORAGE  CELLS 

C 

b 

2 

(COSINECETA))2'5 

Cl 

b 

2 

INTERMEDIATE  VARIABLE 

C2 

b 

2 

INTERMEDIATE  VARIABLE 

CDKAG  170A 

-170C 

1 

SEE  INPUT  LISTING 

CELIP 

A. 30 

INTERMEDIATE  VARIABLE 

COASC 

A. 5b 

INTERMEDIATE  VARIABLE 

COET 

A. 54 

INTERMEDIATE  VARIABLE 

COINCL 

x> 

in 

• 

< 

INTERMEDIATE  VARIABLE 

CON (1-3) 

ASSIGNED  VARIABLES 
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VARIABLE 

EliUA 

REF 

DEFINITION 

C0N(4-b) 

122-139 

1 

INTERMEDIATE  VARIABLES 

CON (6-7) 

17 

1 

INTERMEDIATE  VARIABLES 

CON (8) 

ASSIGNED  VARIABLE 

CON (9-10) 

141  » 
147 

1 

INTERMEDIATE  VARIABLES 

COPER I 

A*  5b 

INTERMEDIATE  VARIABLE 

COSALP 

101 

1 

INTERMEDIATE  VARIABLE 

COSEN 

A.1BA 

INTERMEDIATE  VARIABLE 

LOTAU2 

5 

? 

INTERMEDIATE  VARIABLE 

CSAS1 

98A 1 9flb 

1 

INTERMEDIATE  VAKlARLE 

CT(  9) 

see  Input  listing  (these  are  recomputed  and 

PRINTED  OUT  AT  THE  END  OF  FACH  FILTERING  PHASE 
SEE  EQUATION  A. 160) 

CTAU 

b 

? 

(COSINE ( TAU/2 ) ) 2,5 

LTIN(6) 

A.  16U 

A  •  lf>C 

3.0#4.0*6.0,fl.0#11.0  AND 
.95 

CTNEm (9) 

A. 14 

INTERMEDIATE  VARIABLES 

UATES<3) 

STORAGE  CELLS 

JATSEC(I) 

1-1*2. 

DJUPRT ( I ) *B64U0 • 0  ♦  TIMSEC(I).  TIME  IN  SEC¬ 
ONDS  CORRESPONDING  TO  INITIAL  AND  FINAL 

PRINT  TIMES 

uayflr 

MODIFIED  JULIAN  DATE  OF  THE  OBSERVATION  TIME 
(USED  DURING  THE  EPHEMERIS  COMPUTATION) 

UAYJL 

MODIFIED  JULIAN  DATE  OF  THE  EPHEMERIS  PRINT 
time»output 
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VARIABLE 

EGUA 

REF 

DEFINITION 

DAYREF 

MODIFIED  JULIAN  DATE  OF  RI6HT  ASCENSION  AND 
DECLINATION  OF  SUN 

OAYS 

DAY  OF  YEAR  ON  OBSERVATION  CARD 

DCRATE 

DECLINATION  RATE  (DEG/SEC) »  OUTPUT 

UECLIN 

127 

1 

DECLINATION  (DEG) 'OUTPUT 

DELAT(l) 

1-1 '  • • • '  NWS 

19 

1 

GEODETIC  LATITUDE  -  GEOCENTRIC  LATITUDE  (FOR 
EACH  STATION) 

OELRT 

REFRACTION  CORRECTION  FOR  RANGE 

DELS 

25 

1 

REFRACTION  CORRECTION  FOR  ELEVATION 

OELTAM 

A. 20 

INTERMEDIATE  VARIABLE 

OELTOR (6) 

CORRECTION  TOLERANCES  FOR  X'Y'Z'X' Y'Z. 

THESE  SIX  VARIABLES  ARE  EQUAL  TO  10.0' 10. O' 
10.0' 0.01 ' 0*01 '0.01 'RESPECTIVELY 

OENL ( 3 ) 

ARRAY  OF  INTERMEDIATE  VARIABLES 

OENTEM 

DENSITY  IN  G/CM**3 

UIA6(6»6) 

ASSIGNED  VARIABLES 

DIFR(7) 

49 

1 

DIFFERENCES  BETWEEN  MEASURED  AND  COMPUTED 
OBSERV AT I ONS ' OUTPUT . 

THESE  SEVEN  CELLS  REPRESENT  -  ELEVATION  OR 
DECLINATION' AZIMUTH  OR  RIGHT  ASCENSION 'RANGE 
RANGE  RATE 'ELEVATION  RATE 'AZIMUTH  RATE  AND 
RANGE  ACCELERATION ' RESPECTIVELY 

UJO 

MODIFIED  JULIAN  DATE  OF  JANUARY  l'MOBS(2)  — 
WHERE  MOBS (2)  IS  THE  YEAR  OF  THE  OBSERVATION 

djbgn 

STORAGE  CELL 

□JREF 

CURRENT  EPOCH  OF  THE  BASIC (SIDEREAL)  SYSTEM 

DJREG 

MODIFIED  JULIAN  DATE  CORRESPONDING  TO  THE  RE* 
GULAR  PRINT  TIME 
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VARIABLE 

EQUA 

REF 

DEFINITION 

OJULN 

CURRENT  MODIFIED  JULIAN  DATE 

DJULO 

PREVIOUS  MODIFIED  JULIAN  DATE 

OJUPRT ( 2 ) 

MODIFIED  JULIAN  DATES  OF  INITIAL  AND  FINAL 
PRINT  TIMES 

OLOG 

INTERMEDIATE  VARIABLE 

DNP (37*89) 

TABLE  OF  COMMON  LOGARITHMS  OF  DENSITY  CG/CM3) 
AS  A  FUNCTION  OF  EXOSPHERIC  TEMPERATURE  AND 
ALTITUDE 

UPRINS 

EQUIVALENT  TO  DPR I NT 

OPRINT 

SEE  INPUT  LISTING 

UPRR 

INTERMEDIATE  VARIABLE 

0SAVE(3) 

VALUE  OF  THE  UPDATED  VELOCITY  VECTOR  (X*Y*Z) 

OSPI(I) 

l=l*...*NOSPRI 

MODIFIED  JULIAN  DATES  OF  THE  SPECIAL  PRINT 
TIMES 

OSTEP 

102 

1 

POSITION  PERTURBATION  =.4(KM)#IN  STATE 
TRANSITION  MATRIX 

DTH1 

194 

1 

INTERMEDIATE  VARIABLE 

OUM ( 3 ) 

INTERMEDIATE  VARIABLES 

OV(3) 

POSITION  VECTOR  X*Y*Z  (KM) 

UVEQU2) 

OVEQ( 1 ) »DVEQ(3) »DVEQ(5)  ARE  EQUIVALENT  TO 

DV ( 1 ) * D V  ( 2 ) #  D V  C  3 ) * RESPECT I VEL Y  • 

DVEQ ( 7 ) » DVEQ ( 9 )  • DVEQ (11)  ARE  EQUIVALENT  TO 

DVP ( 1 ) * DVP ( 2 ) * D VP ( 3 ) » RESPECT I VELY 

DVH(3) 

FINAL  INTEGRATED  VALUES  OF  X*Y*Z  IN  INTE¬ 
GRATION  SUBROUTINE 

UVH2(3) 


INTERMEDIATE  VALUES  OF  X*Y*Z  IN  INTEGRATION 
SUBROUTINE 


t 


VARIABLE 

EOUA 

REF 

DEFINITION 

UVH2PC3) 

INTERMEDIATE  VALUES  OF  XtY'Z  IN  INTEGRATION 
SUBROUTINE 

UVHPO) 

FINAL  INTEGRATED  VALUES  OF  X'YtZ  IN  INTE¬ 
GRATION  SUBROUTINE 

DVHPPC3) 

INTERMEDIATE  VALUES  OF  X'Y'Z*  IN  INTEGRATION 
SUBROUTINE 

DV0(3) 

INITIAL  VALUES  OF  X*Y»Z  AT  START  OF  INTE¬ 
GRATION  SUBROUTINE 

UVP(3) 

VELOCITY  VECTOR  X»f»Z  (KM/SEC) 

OVPOO) 

INITIAL  VALUES  OF  X'Y'Z  AT  START  OF  INTE¬ 
GRATION  SUBROUTINE 

UVP00(3) 

INTERMEDIATE  VALUES  OF  XtY'Z  IN  INTEGRATION 
SUBROUTINE 

OVPPO) 

INTERMEDIATE  VALUES  OF  X»Y#Z  IN  INTEGRATION 
SUBROUTINE 

ECA 

91A-97B 

1 

ECCENTRIC  ANOMALY  El  OR  E2 

ECAD 

91A-97B 

1 

ECCENTRIC  ANOMALY  El  OR  E2 

ECANOM 

A. 42 

1 

ECCENTRIC  ANOMALY  USED  DURING  EPHEMERIS 
COMPUTATION 

ECC 

87 

1 

ECCENTRICITYrE 

ECC2 

ECC*ECC 

ECCEN 

SEE  INPUT  LISTING 

ECCENT 

A.31 

ECCENTRICITY » OUTPUT 

eelip 

A. 31 

COMPUTED  ECCENTRICITY  AT  THE  EPOCH  OF  THE 
INPUT  ORBITAL  ELEMENTS 

EELIP2 

EELIP*EEL1P 

EHM 

ALTITUDE  ABOVE  WHICH  DENSITY  IS  COMPUTED  BY 
EXPONENTIAL  EXTRAPOLATION  =  1000*0  KM 
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VARIaBLF 

EoUa 

RFF 

Dfc  FINi TIOn 

tLtMAX 

StF  INPUT  LISTING 

tLbMIN 

sfe  input  listing 

tLtVAT 

1?U 

1 

ELEVATION  U)Eo) # OUTPUT 

U.KAIE 

1^ 

1 

ELEVATION  RATt  ( DtG/SF C )» OUTPUT 

tNDUAT 

RLf)  WORD  -FNDuAT  (StF  INPUT  LISTING) 

t_Pb 

17M 

1 

MEAN  URLIUUITY  OF  FCL1PTIC 

ciRASlX  1  n ) 

ASSIGNED  variables 

cRASt(^O) 

ASSIGNED  variables 

tT 

COMPUTED  tCCENTRlC  ANOMALY  AT  THF  EPOCH  OF 

ThE  Input  orbital  elemfnts 

LTA(o) 

INTERMEDIATE  VARIABLES 

C.TASUN 

1  flrt 

1 

OiRtCllON  COSINE  OF  SUM 

C.TM 

FXOSPMFRIC  TEMPERATURE  BFLOW  WHICH  DENSITY 

computed  oy  extrapolation  =  goo.o  peg  k 

IS 

tXTRFM 

INCREMENT  of  LXOSPHER1C  TtMpERATURF  BELOW 
LOWEST  TEMPERATURE  OF  DNP  DATA  ALLOWED  FOR 
temperature  branch  =  ino.u  deo  k 

LO 

F 

o 

1 

StF  INPUT  LISJING 

H(b) 

lSrt-1 SH 

1 

intermediate  variables 

F  C  A 

1 

intermediate  variable 

hCUMST 

170A- 

l7i»C 

1 

intermediate  variable 

i-Earih 

n 

1 

(  1  .  U  -  F  )2>  Fe 

F ILHlRl?) 

bld  information  used  for  hfadimg  of  each 
page 
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VARIABLE 

EQUA  REF  DEFINITION 

FLO 

0-0 

FL1 

1.0 

FL10 

10*0 

FL12 

12*0 

FL1P5 

1.5 

FL2 

2*0 

FL24 

24.0 

FL2PI 

2.0*PI 

FL3 

3.0 

FL3600 

3600. Or NUMBER  i 

FL3652 

36525 .Or  NUMBER 

FL4 

4.0 

FL5 

5.0 

FL6 

6.0 

FL60 

60. 0» NUMBER  OF 

FL7 

7.0 

FL8 

8.0 

FL864H 

86400 • 0  *  NUMBER 

FL9 

9.0 

FL96 

96.0 

FLP001 

0.001 

FLP01 

0.01 

FLP1 

0.1 
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VARIABLE 

FnIIA 

REF 

ntFINiTIOu 

HLPS 

O.S 

ELPI 

Pi  #3.14l5y?6*S 

RLRAlj 

.til74S3?92S*CoNVEl<SI0N  FACTOR  OEGRFES  TO 
RADIANS 

RSM1 

IS« 

1 

INTERMEDIATE  vARImBLF 

ESM2 

isr 

1 

INTERMEDIATE  vAkIaRLF 

o 

4 

2 

.78S30A16  RADIANS  (UNUSED) 

oAm 

1  7rt 

1 

LONG 11 DUE  OF  bUN 

oAmS 

24 

1 

TOTAL  REFRACTION  BENDING  1HR0UGH  THE  TRO- 
POSHFkF  (hAOImNS) 

oEOCEN 

A  •  S  / 

GEOCEnTrIC  laiituuf  (UFG) »outpht 

oEODtT 

A. 61 

GEODETIC  LATIIUUF  (UFG) » OUTPUT 

M(NNTP#6) 

122- 

l4fl 

1 

H  MAT»<  I  x  #  THE  MATRIX  OF  THt  PARTIAL  DERIVA¬ 
TIVES  WITH  RE SPEC  r  TO  THE  SIX  ORBITAL 

FlFMEnTS  IN  TmE  FuRM  OF  POSITION  AND  VELOCITY 

coordinates 

mCuEFF 

4 

2 

intermediate  variable 

mCOH  (6) 

23 

1 

increment  laylrs  used  in  rffraCtiou  Cor¬ 
rect  ion  SUBROUTINE 

riUBAb 

4 

2 

DaFFEhEnCl  OF  RaSaT  -  SUNK A 

Ml  NO 

102 

1 

VARIABLE  INTEGRATION  INTERVAL ( SEC . ) 

MIND2 

(hi  NO  )* 

MlNDiiB 

2 

(hi NO ) /B • U 

MlNO^H 

(hII\l()//B,0 
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VARIABLE 

EQUA 

REF 

DEFINITION 

hlND29 

(HlND)2/96.0 

HIND02 

H1N0/2.0 

HI NO 06 

HIN0/6.0 

HIN024 

H1NU/24.0 

HNOW 

ACCUMULATIVE  INTEGRATION  INTERVAL  (THE 

LIMIT  IS  HTOTAL) 

HPHTR(7#7) 

ASSIGNED  VARIABLES.  MATRIX  TO  BE  INVERTED 

OR  ITS  INVERSE. 

rtRPRTU) 

SEE  INPUT  LISTING 

HTOFLR 

INTEGRATION  INTERVAL  BETWEEN  TWO  OBSERVATION 
TIMES (SEC.) 

HT0T1 

9b 

1 

INTERMEDIATE  VARIABLE 

HTOTAL 

INTEGRATION  INTERVAL  BETWEEN  TWO  PRINT  TIMES 
OR  BETWEEN  TWO  OBSERVATION  TIMES(SEC.) 

HTOTAS 

EQUIVALENT  TO  HTOTAL 

HTR<6»NNTP) 

5U 

1 

THE  WEIGHTING  MATRIX  OR 

50 

1 

THE  TRANSPOSE  OF  H 

1 

PROGRAM  INDEX 

112345 

PROGRAM  LOGIC  CONTROL 

1AM 

INTERMEDIATE  VARIABLE 

AB<7) 


STORAGE  CELLS 


e 


VARIABLE  EQUA  REF  DEFINITION 

1COUNT  ASSIGNED  VARIABLE 

ICT  ITERATION  COUNTER 

IDIR  =1*  FILTER  IN  BACKWARD  DIRECTION 

=0*  FILTER  IN  FORWARD  DIRECTION 

i equal  program  logic  control 

XEX  EXPONENT  OF  RANGE  ON  AN  OBSERVATION  CARD 

IFF (7)  PROGRAM  INDICES  AND  STORAGE  CELLS 

IFILSV  PREVIOUS  VALUE  OF  IFILTR 


IFILTR  =1*  DJULN  AND  TNORMN  CORRESPOND  TO  AN  EPHE- 

MER1S  PRINT  TIME 

=0’  DJULN  AND  TNORMN  CORRESPOND  TO  AN  OBS¬ 
ERVATION  TIME 

=-l*DJULN  AND  TNORMN  CORRESPOND  TO  AN  EPHE- 
MERIS  AND  AN  OBSERVATION  TIME 

1FINIS  =1 'FINAL  EPHEMERIS  PRINT  OUT 

=0*EPHEMERIS  PRINT  OUT  TO  CONTINUE 

IFIRST  =0*1  FILTER  ONLY 

=2*  FILTER  AND  PRINT  EPHEMERIS 

1G(7)  PROGRAM  INDICES  AND  STORAGE  CELLS 

AGO  =-l 'TRANSFORMATION  OF  COVARIANCE  MATRIX 

RESULTED  IN  NEGATIVE  DIAGONAL  ELEMENTS. 
PROGRAM  WILL  NOT  CONTINUE. 

=0*1  TRANSFORMATION  WAS  SUCCESSFUL 
(SEE  EQUATION  51*  REFERENCE  1) 

II  PROGRAM  INDEX 


I IN  SYSTEM  INPUT  TAPE 

1J  NUMBER  OF  STATIONS  IN  A  BINARY  DATA  RECORD 

OBSERVING  THE  SATELLITE ( IJ *  10) 
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.4:  i  AiAaMMWM? 


I 


VARIABLE  EQUA  REF  DEFINITION 

1JFILP  =1  • INCREMENT  JFILP  BY  1 , I. E. FILTER  FORWARD 

=-l* INCREMENT  JFILP  BY  -1 * I •£. FILTER  BACKWARD 

1JKL  STORAGE  CELL*  EQUALS  IOUT  OR  MBCD 

(SEE  IOUT  AND  MBCD) 

1L  =1*EPHEMERIS  PRINT-OUT  IS  AT  A  SPECIAL 

PRINT  TIME 

=0*EPHEMERlS  PRINT-OUT  IS  AT  A  REGULAR 
PRINT  TIME 

ILIMIT  =0*DIFR  ARI.  WITHIN  THE  TOLERANCES 

=1*DIFR  ARE  NOT  WITHIN  THE  TOLERANCES 
(SEE  DIFR) 

1LMTO  VALUE  OF  ILIMIT  AT  THE  FINAL  OBSERVATION  CARD 

1NCHM  INTEGRAL  VALUE  OF  XINCHM  (INTEGER) 

iNCTM  INTEGRAL  VALUE  OF  XINCTM  (INTEGER) 

INTEND  =-l* INTEGRATION  IS  COMPLETE 

=0*  INTEGRATION  IS  NOT  COMPLETE 

iONE  =0* REWIND  TAPE  MBCD 

=1*00  NOT  REWIND  TAPE  MBCD 

IOUT  SYSTEM  OUTPUT  TAPE 

1PAGE  PAGE  COUNTER  FOR  FILTERING  OUTPUT 

I  PHI  PROGRAM  INDEX 


1REC  PROGRAM  INDEX 

iREV  REVOLUTION  NUMBER  *  OUTPUT 

1REV0  SEE  INPUT  LISTING 

1REWD  =0*REWIND  TAPE  MBIN 

=1*00  NOT  REWIND  TAPE  MBIN 

1SAVE  PROGRAM  INDEX 
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VARIABLE 


EtiUA  REF  DEFINITION 


1SCLH 

ISH 

iSMOOH 

1ST 

1TLT 

iTIMt 

1TYPE 


J 

JACK 

JB 

JDAPKT12) 

JERR 

JF1LP 

JFiRbT 


Jj 

JM 

JMOPRT (2) 


INTEGRAL  value  OF  XISCLH  /  XINCHM  (INTEGER) 
INTEGRAL  VALUE  OF  SHM  (INTEGER) 

SEE  INPUT  LISTING 

INTtGKAL  VALUE  OF  STM  (INTEGER) 

ITERATION  COUNTER 
STORAGE  CELL 

=1  *P INAL  DATA  RECORD  AT  AN  EPHEMERIS 
PRINT  TIME 

=-l» ANOTHER  DATA  RECORD  TO  FOLLOW  AT  AN 
EPHEMERIS  PRINT  TIME 

program  index 

=U. STORE  PROGRAM  CONSTANTS 
=1 .CALCULATE  INITIAL  VALUES 

PROGRAM  INDEX 

SEE  INPUT  LISTING 

PROGRAM  LOGIC  CONTROL 

OBSERVATION  CARD  COUNTER 

=0 .READ  OBSERVATION  CARDS  IN  INCREASING  ORDER 
ANu  FILTER 

=1 .READ  OBSERVATION  CARDS  IN  DECREASING  ORDER 
ANU  FILTER 

=2. READ  OBSERVATION  CARDS  IN  INCREASING  ORDER 
.FILTER  AND  PRINT  EPHEMERIS 

PROGRAM  INDEX 

program  index 
see  input  listing 


VARIABLE 


EttUA  REF  DEFINITION 


JSP 

JSPDA 

uSPMO 

JSPSAV 

JSPYR 

JTYPRT 

JYRPRK2) 

K 

KAR6 

KARGSV 

KBEG 

KCLASS 

KCOUNT 

KDaOUT 

KOECOB 

KEY 

KEYTAP 

KFO 

KF1 

KF10 

KF100 

KF2 


PROGRAM  COUNTER  FOR  SPECIAL  PRINT  TIMES 
StE  INPUT  LISTING 
StE  INPUT  LISTING 

VALUE  OF  JSP  AT  THE  PREVIOUS  OBSERVATION  TIME 

StE  INPUT  LISTING 

SEE  INPUT  LISTING 

SEE  INPUT  LISTING 

PROGRAM  INDEX 

PROGRAM  LOGIC  CONTROL 

VALUE  OF  KARG  AT  THE  OBSERVATION  TIME 

PROGRAM  INDEX 

SATELLITE  CLASSIFICATION  ON  OBSERVATION  CARD 
StE  INPUT  LISTING 
CALENDAR  DAY# OUTPUT 

DECADE  OF  TIME  IN  YEARS C50 »60»70» •  •• ) 

PROGRAM  INDEX 
StE  INPUT  LISTING 
0 
1 

10 

100 

2 
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T 


variable 

KF25 

KF  3 

KF3i 

KF4 

KFb 

KFb 

*FbO 

KF7 

KFtt 

KF9 

KFILPR 

KFILPT 

KF1N 

KHROUT 

kiter 


klamra 

KLONUS 

KMIOUT 

KMOOUT 


EQUA  REF  DEFINITION 
2b 

3 
31 

4 

5 

6 
60 

7 

8 
9 

=0 * PROGRAM  IS  FILTERING  AND  COMPUTING 
EPHEMER IS 

si 'PROGRAM  IS  COMPUTING  ONLY  THE  EPHEMERIS 

so 'PROGRAM  IS  FILTERING 

si 'PROGRAM  IS  FILTERING  AND  COMPUTING 

EPHEMERIS 

PROGRAM  INDEX 

HOUR  OF  DAY 'OUTPUT 

si 'OBSERVATION  CARD  HAS  BEEN  SAVED  FOR 
SMOOTHING 

so 'OBSERVATION  CARD  HAS  NOT  BEEN  SAVED  FOR 
SMOOTHING 

PROGRAM  LOGIC  CONTROL 
PROGRAM  INDEX 
MINUTE  OF  HOUR 'OUTPUT 
CALENDAR  MONTH' OUTPUT 
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variable 


EQUA  REF  DEFINITION 


KOBSPR 

KONCE 


KQ 

*Q20 

MOT 

MOTM1 

MYPRT 

KYROUT 

L(7) 

LASTCA 

LFO 

LINES 

LINOUT 

LMBDA 

LOOKAN 

LPAGE 

LUP 

LUV 

LUV1 

MBCD 


SEE  INPUT  LISTING 

=0*THE  FIRST  SIX  COR  LESS)  OBSERVATION  CARDS 
HAVE  BEEN  SAVED  FOR  SMOOTHING 
=1* OBSERVATION  CARDS* OTHER  THAN  THE  FIRST  SIX 
HAVE  BEEN  SAVED  FOR  SMOOTHING 

PROGRAM  INDEX 

INTERMEDIATE  VARIABLE 

PROGRAM  INDEX 

PROGRAM  INDEX 

(SEE  JTYPRT) 

CALENDAR  YEAR* OUTPUT  CLAST  2  OIGITS  OF  19XX) 

PROGRAM  INDICES  AND  STORAGE  CELLS 

LAST  VALUE  OF  JFILP  AT  WHICH  AN  OBSERVATION 
CARD  HAS  BEEN  ACCEPTED 

LINE  COUNTER  INCREMENT 

LINE  COUNTER  FOR  FILTERING  OUTPUT 

LINE  tOUNTER  FOR  EPHEMERIS  OUTPUT 

PROGRAM  INDEX 

PROGRAM  LOGIC  CONTROL  FOR  EPHEMERIS  OUTPUT 
PAGE  COUNTER  FOR  EPHEMERIS  OUTPUT 
PROGRAM  INDEX 
PROGRAM  INDEX 
PROGRAM  INDEX 

BCD  OUTPUT  TAPE  REQUESTED  BY  CUSTOMER 
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I 


VARIABLE 

EQUA 

REF 

DEFINITION 

MBXN 

BINARY  OUTPUT  TAPE  REQUESTED  bY  CUSTOMER 

ME(J(  J) 

U-l  ,  •  •  •  9  9 

=U » COMPUTE  TRANSFORMATION  MATRIX  TO  THE  MEAN 
EQUINOX  ANU  EQUATOR  OF  THE  CELESTIAL  REF¬ 
ERENCE  SYSTEM 

=1# TRANSFORMAT I ON  MATRIX  HAS  BEEN  COMPUTED 
(SEE  PPRIME ( 3,3, J) 

MEKASE12U) 

ASSIGNED  VARIABLES 

MMTAPE 

STOhAgE  TAPE  USED  BY  PROGRAM 

MOBS ( 1 ) 

SATELLITE  NUMBER  ON  OBSERVATION  CARD 

MOBS (2-4) 

SEE  INPUT  LISTING 

OR  YEAR, MONTH, DAY  OF  OBSERVATION 

MOLO 

STORAGE  CELL 

MSTHR 

HOUR  OF  GREENWICH  MEAN  SIDEREAL  TIMErOUTPUT 

MSTMIN 

MINUTE  OF  GREENWICH  MEAN  SIDEHEAL  TIME, 
OUTPUT 

N 

PROGRAM  INDEX 

N2 

A.23C 

integral  value,  of  bn 

A  .230 

INTEGRAL  VALUE  OF  BNO 

NAT 

SEE  INPUT  LISTING 

NCARSV 

TOTAL  NUMbER  OF  OBSERVATION  CARDS 

NCOL 

COLUMN  ERROR  DESIGNATOR 

ND1 

20 

1 

=2  »  THE  MATRIX  PITH (SEE  P1TR)  HAS  BEEN  RECOM¬ 
PUTED,  HENCE  PN1  MUST  BE  RECOMPUTED 
=1»P1TR  MATRIX  HAS  NOT  BEEN  RECOMPUTED 

NDIM 

FORTRAN  DIMENSION  OF  A  SQUARE  MATRIX 
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VARIABLE 

NDIR 

NE2 

NEQ 

NE01 

NF 

NFIRST 

NG 

NHM 

NHMSCL 

NMS 

NNEQ 

NNTP 


NO 

NOB 

NOSAT 

NOSPRI 


EGUA  REF  DEFINITION 

=1# FILTER  AND  COMPUTE  EPHEMER1S 
=2 ♦ COMPUTE  EPHEMERlS»NO  FILTERING 
=3* COMPUTE  EPHEMERIS'NO  FILTERING .SET 
ND1R=1 

YEAR  OF  OBSERVATION 

CODE  OF  THE  CELESTIAL  REFERENCE  SYSTEM  ON 
THE  OBSERVATION  CARD 

STORAGE  CELL  FOR  NEQ 

PROGRAM  INDEX 

INITIALIZATION  SECTION  CONTROL  (INTEGER) 
PROGRAM  INDEX 

NO.  OF  ALTITUDES  IN  DNP  TABLE  =  89  (INTEGER) 

INDEX  FOR  HIGH  ALTITUDE  BRANCH 

SEE  INPUT  LISTING 

PROGRAM  INDEX 

IF  NTP  =1#NNTP  =2 
IF  NTP  =1*  NNTP  =NTP 

NNTP  REPRESENTS  THE  NUMBER  OF  SIMULTANEOUS 
OBSERVATIONS  ON  AN  OBSERVATION  CARD. 

NNTP  IS  ALSO  USED  AS  THE  DIMENSION  OF  A  MATRIX 
TO  BE  INVERTED (SEE  HPHTR ) 

STORAGE  CELL 

SEE  INPUT  LISTING 

SEE  INPUT  LISTING 

SEE  INPUT  LISTING.  NOSPRI  IS  DECREASED  BY  1 
EACH  TIME  THE  PROGRAM  PRINTS  AT  A  SPECIAL 
PRINT  TIME 
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« 


variable 


EQUA  REF  DEFINITION 


NOSPRS 

NRt 

NS 

NSG(l) 

i=l*...*NMS 

NSTPHT 

NTIME 

NTM 

NTP 


NTPOLD 

NUMBER (111) 

NUMSTA ( I ) 
1=1*...*NMS 

NWORQS 

OBSC(7) 

0BSM(7) 


VALUE  OF  NOSPNl  AT  THE  PREVIOUS  OBSERVATION 
TIME 

=1* INITIALIZE  PROGRAM  LOGIC  FOR  COMPUTATION 
OF  DECLINATION  AND  RIGHT  ASCENSION 
=2* INITIALIZATION  HAS  BEEN  DONE 

PROGRAM  INDEX 

16*17  1  SEE  INPUT  LISTING 


SEE  INPUT  LISTING 
STORAGE  CELL 

NO.  OF  TEMPERATURES  IN  DNP  TABLE  =  37  (INTEGER) 

CODE  WHICH  DESIGNATES  TYPE  OF  SIMULTANEOUS 
OBSERVATIONS 

=1  * DECLINATION* RIGHT  ASCENSION 
=2  ELEVATION* AZIMUTH 
r 3  *  ELEVATION  *  AZ I MUTH *  RANGE 
=4  *  ELE V AT I ON  *  AZ I MUTH • R  ANGE  *  RANGE  RATE 
27* SAME  AS  =4  AND  ELEVATION  RATE* AZIMUTH 
RATE*  RANGE  ACCELERATION 

STORAGE  CELL* VALUE  OF  NTP 

EQUIVALENT  TO  ANGDAT ( 111 ) 

SEE  INPUT  LISTING 


PROGRAM  INDEX 

120-  COMPUTED  OBSERVATIONS. THESE  SEVEN  VARIABLES 

141  REFER  TO  ELEVATION  OR  DECLINATION* AZIMUTH  OR 

RIGHT  ASCENS I ON  t RANGE  *  RANGE  RATE* ELEVATION 
RATE* AZIMUTH  RATE  AND  RANGE  ACCELERATION  *  RES¬ 
PECTIVELY  • 

MEASURED (INPUT)  OBSERVATIONS  ON  THE  OBSER¬ 
VATION  CARD. THESE  SEVEN  VARIABLES  ARE  ANALO¬ 
GOUS  TO  ObSC ( 7) • 


VARIABLE 

EQUA 

REF 

definition 

OBSNO 

COLUMNS  73-76  OF  OBSERVATION  CARD 

OINCL 

SEE  INPUT  LISTING 

OINCLI 

A. 37 

INCLINATION  (DEG) 'OUTPUT 

OKINT1 

193 

1 

INTERMEDIATE  VARIABLE 

OLAMO 

A. 64 

SATELLITE  LONGITUDE  (DEG) 'OUTPUT 

OLAMS 

A.64 

SATELLITE  LONGITUDE  (RAD) 

OLAMSP 

A. 64 

SATELLITE  LONGITUDE  (HAD) 

OLODAY 

VALUE  OF  MODIFIED  JULIAN  DATE (DJULN)  AT  WHICH 
AN  OBSERVATION  CARD  HAS  BEEN  ACCEPTED 

OLDTIM 

VALUE  OF  TNORMN  CORRESPONDING  TO 

DJULN  (SEE  OLDDAY) 

omeana 

A. 43 

MEAN  ANOMALY  (DEG) 'OUTPUT 

OMS 

180 

1 

INTERMEDIATE  VARIABLE 

omsat 

StE  INPUT  LISTING 

OMU 

SEE  INPUT  LISTINGS 

ONE 

1.0 

ONUM 

INTERMEDIATE  VARIABLE 

P 

4 

2 

•20943951  RADIANS 

PO 

INTERMEDIATE  VARIABLE 

Pl2(3*3) 

ASSIGNED  VARIABLES 

P1TR(3'3) 

15*20 

1 

transpose  of  nutation-precession  matrix 

P2(3»3) 

ASSIGNED  VARIABLES 

P2TR(3'3) 


lb  1 


transpose  of  nutation-precession  matrix 

ASSIGNED  VARIABLES 

NUTATION-PRECESSION  MATRIX  (COMPUTED  AT  DJULN) 


variable 

EOUA 

REF 

DEFINITION 

HASS 

SEE  INPUT  LISTING 

PERlGE 

A.  41 

ARGUMENT  OF  PERIGE  (DEG). OUTPUT 

PERIOD 

PERIOD  OF  SATELLITE  AT  THE  EPOCH  OF  THE 

INPUT  ORBITAL  ELEMENTS (SEC ) 

PERMUT (7) 

STORAGE  CELLS 

PHI  loro)  51 

*106 

1 

STATE  TRANSITION  MATRIX 

PH1I (7*7) 

5U 

1 

INTERMEDIATE  MATRIX 

PH1LAT ( I ) 
1=1*...*NMS 

SEE  INPUT  LISTING 

PH1TR ( o»b) 

ASSIGNED  VARIABLES 

PMAT \6*6) 

50- 

52 

1 

COVARIANCE  MATRIX 

PN1 (3*3) 

20 

1 

TRANSFORMATION  MATRIX  FROM  THE  TRUE  EQUINOX 

and  equator  of  date  to  the  mean  equinox  and 

EQUATOR  OF  THE  PARTICULAR  CELESTIAL  SYSTEM 

PNlTR (4*3) 

129* 

131 

1 

TRANSPOSE  OF  THE  PN1  MATRIX 

PNODAL 

A. 25 

NODAL  PER 100  OF  SATELLITE  AT  THE  EPOCH 

OF  THE  INPUT  ORBITAL  ELEMENTS  (SEC.) 

PPRIME(3*3* J) 
J-l * i *  *  *9 

20 

1 

TRANSFORMATION  MATRICES  TO  THE  MEAN  EQUINOX 

AND  EQUATOR  OF  THE  CELESTIAL  REFERENCE  SYSTEM 

U(7) 

55 

1 

StE  INPUT  LISTING  (THESE  ARE  RECOMPUTED  AND 
PRINTED  OUT  AT  THE  END  OF  EACH  FILTERING  PHASE 
SEE  EQUATION  A. 15) 

OP  (2) 

55 

1 

SEE  INPUT  LISTING  (THESE  ARE  RECOMPUTED  AND 
PRINTED  OUT  AT  THE  END  OF  EACH  FILTERING  PHASE 

see  Equation  a.is) 

H 

5 

2 

•  3 
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VARIABLE 

EGUA 

REF 

DEFINITION 

R2 

99 

1 

INTERMEDIATE  VARIABLE 

RADIUS 

DISTANCE  OF  THE  SATELLITE  FROM  THE  CENTER  OF 
THE  EARTH  AT  THE  EPOCH  OF  THE  INPUT  ORBITAL 
ELEMENTS 

RADPRT 

A. 24 

DISTANCE  OF  THE  SATELLITE  FROM  THE  CENTER  OF 
THE  EARTH  (KM)*OUTPUT 

RALON 

18 

1 

RIGHT  ASCENSION  OF  THE  OBSERVING  STATION 

RANGES 

140 

1 

RANGE  (KM) * OUTPUT 

RANRAT 

141 

1 

RANGE  RATE  (KM)* OUTPUT 

RARATE 

RIGHT  ASCENSION  RATE  (DEG/SEC) * OUTPUT 

RASAT 

RIGHT  ASCENSION  OF  SATELLITE  (RADIANS) 

RASE (6) 

2 

CORRECTION  IN  POSITION  AND  VELOCITY  VECTOR 

AT  THE  LAST  OBSERVATION  CARD  WHICH  HAS  BEEN 
ACCEPTEO. THESE  ARE  COMPUTED  WHENEVER  THE 
PROGRAM  ITERATES* OUTPUT 

RDPP 

170A- 

170C 

1 

INTERMEDIATE  VARIABLE 

REARTH 

7 

1 

SEE  INPUT  LISTING*RC 

RER 

158* 

159 

1 

INTERMEDIATE  VARIABLE 

HER  2 

158* 

159 

1 

INTERMEDIATE  VARIABLE 

RESULT 

INTERMEDIATE  VARIABLE 

RFC (6) 

INTERMEDIATE  VARIABLES 

RHI (7*7) 

50* 

52 

1 

INTERMEDIATE  MATRIX 

RHO  170A 

-170C 

1 

ATMOSPHERIC  DENSITY  (KG/M3) 
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T 


VARIABLE 

EUUA 

REF 

DEFINITION 

RMEAN 

REARTh*  ( 2  •  0-F )  /2  •  0 

KT2 

12b 

1 

INTERMEDIATE  VARIABLE 

KTASC 

12b 

1 

RIGHT  ASCENSION  <DEG)*OUTPUT 

HTC2 

122*134- 

147 

1 

RANGE  SQUARED 

HX 

128 

1 

INTERMEDIATE  VARIABLE 

KX2 

129 

1 

INTERMEDIATE  VARIABLE 

S 

b 

2 

(SINE C THETA)  )*‘S 

bAl 

bOA 

1 

INTERMEDIATE  VARIABLE 

bA2 

bOd 

1 

INTERMEDIATE  VARIABLE 

SA3 

80L 

1 

INTERMEDIATE  VARIABLE 

SAPA 

88A 

1 

INTERMEDIATE  VARIABLE 

bATRAD 

DISTANCE  OF  THE  SATELLITE  FROM  THE  CENTER  OF 
THE  EARTH  (KM) 

bAvE (3) 

VALUE  OF  THE  UPDATED  POSITION  VECTOR  (X*Y*Z) 

bD 

80A- 

80C 

1 

INTERMEDIATE  VARIABLE 

bDV (b*b) 

VALUES  OF  POSITION  AND  VELOCITY  VECTOR  AT 
six  observation  TIMES.  THESE  MILL  BE  USED  FOI 
SMOOTHING. 

bECOUT 

SECONDS  OF  MINUTE* OUTPUT 

bECPMT (2) 

StE  INPUT  LISTING 

bHM 

LOWEST  ALTITUDE  IN  DNP  TABLE  =  120.0  KM 

blASC 

A. 5b 

INTERMEDIATE  VARIABLE 

55 


VARIABLE 

EQUA 

REF 

definition 

SIET 

A  *54 

INTERMEDIATE  VARIABLE 

SIGNA 

=1.0 » VALUE  OF  AN  OBSERVATION  INPUT  IS  POSITIVE 
=-1*0* VALUE  OF  AN  OBSERVATION  INPUT  IS  NEGATIVE 

SlGXX(b) 

A.lA- 

A.2C 

ERRORS  IN  POSITION  ANO  VELOCITY  DUE  TO  A 
STANDARD  timing  ERROR 

SIINCL 

A. 5b 

INTERMEDIATE  VARIABLE 

SINALP 

101 

1 

INTERMEDIATE  VARIABLE 

SINE (6) 

26 

1 

INTERMEDIATE  VARIABLES 

SINEN 

A  *188 

INTERMEDIATE  VARIABLE 

SIPA 

88b 

1 

INTERMEDIATE  VARIABLE 

blPERl 

A. 5b 

INTERMEDIATE  VARIABLE 

bLAT ( I ) 
i=l» • . . »NMS 

b 

1 

GEOCENTRIC  LATITUDE  FOR  EACH  STATION 

bLON ( I ) 

1-1  *  •  •  •  *  NMS 

18 

1 

SEE  INPUT  LISTING 

SMALLA 

A  *32 

COMPUTED  SEMI-MAJOR  AXIS  (KM) 

SMATS(b»b) 

VALUE  OF  PMAT(6'6)  AT  AN  OBSERVATION  WHICH 

HAS  BEEN  ACCEPTED 

SMDAY(b) 

VALUES  OF  DJULN  CORRESPONDING  TO  THE  STORED 
VALUES  OF  THE  POSITION-VELOCITY  VECTOR  (SEE 
SDV(6»6>> 

SMTlM(o) 

VALUES  OF  TNOKMN  CORRESPONDING  TO  SMDAY(6) 

bOA 

84 

1 

INTERMEDIATE  VARIABLE 

bOB 

8b 

1 

INTERMEDIATE  VARIABLE 

bOC 

8b 

1 

INTERMEDIATE  VARIABLE 

56 


VARIABLE 

EGUA 

REF 

DEFINITION 

bOGA 

A. 58 

INTERMEDIATE 

VARIABLE 

bOGB 

A. 61 

INTERMEDIATE 

VARIABLE 

bOGBO 

A. 61 

INTERMEDIATE 

VARIABLE 

SOGC 

A. 60 

INTERMEDIATE 

VARIABLE 

bOK 

90 

1 

INTERMEDIATE 

VARIABLE 

bOKB 

90 

1 

INTERMEDIATE 

VARIABLE 

bOKC 

93- 

1 

INTERMEDIATE 

VARIABLE 

96 

bOKV 

10UA- 

100B 

1 

INTERMEDIATE 

VARIABLE 

bOb 

84 

1 

INTERMEDIATE 

VARIABLE 

bPHR 

StE  INPUT  LISTING 

bPMAT(b»6> 

DIAGONAL  ELEMENTS  OF  THE  COVARIANCE  MATRIX 
(PMAT)  CORRESPONDING  TO  SDV  (SEE  bDV) 

bPMI 

SEE  INPUT  LISTING 

bPSEC 

SEE  INPUT  LISTING 

bQTlME 

A. 54 

INTERMEDIATE 

VARIABLE 

SQTMUA 

A  *54 

INTERMEDIATE 

VARIABLE 

b0T0LR(7) 

53 

1 

SQUARE  ROOT  OF  DIAGONAL  ELEMENTS  IN  R  MATRIX 
THESE  SEVEN  VARIABLES  REFER  TO  ELEVATION 

OR  DECLINATION*  AZIMUTH  OR  RIGHT  ASCENSION' 
RANGE 'RANGE  RATE 'ELEVATION  RATE' AZIMUTH  RATE 
OR  RANGE  ACCELERATION  *  RESPECTIVELY 

bRAO(I) 

7 

1 

GEOCENTRIC  RADIUS  OF  EACH  STATION 

!—!'•••' NMb 


bSAT 


SEE  INPUT  LISTING 


VARIABLE 

EQUA 

REF 

DEFINITION 

»TEP2 

101 

1 

DSTEP  ♦  DSTEP  OR  VSTEP  +  VSTEP 

>TM 

HIGHEST  TEMPERATURE  IN  DNP  TABLE  =  2400*0  DEG 

sTOME 

97A 

1 

INTERMEDIATE  VARIABLE 

bTOU 

INTERMEDIATE  VARIABLE 

aTOUSQ 

97b 

1 

INTERMEDIATE  VARIABLE 

STSECO 

SECONDS  OF  GREENWICH  Ml AN  SIDEREAL  TIME'OUTPUT 

SUM  C  9 ) 

A.  16b 

TOTAL  NUMBER  OF  EACH  TYPE  OF  OBSERVATIONS. 
THESE  REFER  TO  DECLINATION 'RIGHT  ASCENSION' 
ELEVAT ION*  AZIMUTH*  RANGE ' RANGE  RATE  *  ELEVATION 
RATE'AZIMUTH  RATE  AND  RANGE  ACCELERATION 

bUNDEC 

DECLINATION  OF  SUN 

bUNRAS 

RIGHT  ASCENSION  OF  THE  SUN 

I 

S 

2 

EXOSPHERIC  TEMPERATURE  (DEG  K) 

TAbLt 

INTERMEDIATE  VARIABLE 

rc 

13 

1 

MODIFIED  JULIAN  DATE  DIVIDED  BY  36525.0 

TEMl 

DECLINATION  OK  ELEVATION  (DEG)  ON  OBSERVATION 
CARD 

TEMP 

ASSIGNED  VARIABLE 

TEMPt 

ASSIGNED  VARIABLE 

1EMTHD 

90 

1 

INTERMEDIATE  VARIABLE 

TEMTHN 

90 

1 

INTERMEDIATE  VARIABLE 

THtT 

90 

1 

INTERMEDIATE  VARIABLE 

T 1 

TC  (SEE  TC) 

T 12 

TC*TC  (SEE  TC) 

113 

TC*TC*TC  (SEE  TC) 

58 


1 


VARIABLE 

UEQUI(I) 
1=1  * . *  .  *9 


TIMFLR 

1 IMREG 
T I MSEC (2) 

TlMSSSU) 

TIMTOL 

TlNl 

TITLE(il) 

TN2PHI 

TNEXOS 

TNO 

TNORMN 

TNORMO 


EuUA  REF  DEFINITION 

YEARS  OF  STANDARD  CELESTIAL  REFERENCE  SYSTEMS. 

TiEQUim  -  year  of  date  <on  the  observation 
CARD) 

T1EUU1 (2)  -  1900.0 
T1E0UK3)  -  1920.0 
T1EQUK4)  -  1975.0 
T1EGUI (5)  -  2000.0 
TIEGUKb)  -  1650.0 
T1EGU1 (7)  -  1655.0 
T1EQU1 (6)  -  1675.0 
T1EQU1 (9)  -  1960.0 

THE  JULIAN  DAYS  CORRESPONDING  TO  THESE  EPOCHS 
HAVE  bEEN  CONVERTED  TO  MODIFIED  JULIAN  DAYS 
AND  DIVIDED  BY  36525.0 

TIME  OF  DAY  IN  SECONDS  CORRESPONDING  TO 
DAYFLR 

TIME  OF  DAY  IN  SECONDS  CORRESPONDING  TO  DJREG 

time  of  day  in  seconds  corresponding  to  the 

INITIAL  AND  FINAL  PRINT  TIMES (SEE  DJUPRT) 

equivalent  to  timsecu) 

SEE  INPUT  LISTING 

A. 22  D JULN*864U0 • 0  +  TNORMN.  TIME  IN  SECONDS 

CORRESPONDING  TO  THE  INPUT  DATE  AND  TIME 

SEE  INPUT  LISTING 

A. 60  INTERMEDIATE  VARIABLE 

SEE  INPUT  LISTING 

A. 21  TIME (SEC)  TO  THE  FIRST  EQUATORIAL  CROSSING 

FROM  THE  TIME  OK  THE  INPUT  ORBITAL  ELEMENTS 

TIME  OF  DAY  In  SECONDS  CORRESPONDING  TO  DJULN 

TIME  OF  DAY  IN  SECONDS  CORRESPONDING  TO  DJULO 


r 

5 

■ 

V 


VARIABLE 

EoUA 

REF 

definition 

TNORNS 

EQUIVALENT  TO  TNOKMN 

TNOROS 

EQUIVALENT  TO  TNORMO 

1 OLER (7) 

REJECTION  TOLERENCES  (DEG) * OUTPUT 

TOTSEO 

OLDDAY  *  56400.0  +  OLUTIM.  TIME  IN  SECONDS 

Of-  OBSERVATION  LASTCA  (SEE  LASTCA)  WHEN  THE 
PROGRAM  CALLS  SUBROUTINE  ORBINT  FOR  THE  FIRST 
TIME. 

TRIGE(b) 

24-25 

1 

INTERMEDIATE  VARIABLES 

TSEC 

STORAGE  CELL 

rspid) 

1—1 »  .  .  . rNOSPRl 

TIME  OF  DAY  IN  SECONDS  CORRESPONDING  TO  DSPI 

U 

INTERMEDIATE  VARIABLE 

UMlANI 

92 

1 

INTERMEDIATE  VARIABLE*MEAN  ANOMALY 

UMEAN2 

94 

1 

INTERMEDIATE  VARIABLE* MEAN  ANOMALY 

VN 

A. 17 

INTERMEDIATE  VARIABLE 

VSQD 

A. 2b 

INTERMEDIATE  VARIABLE 

VSTEP 

102 

1 

VELOCITY  PERTURBATION  =.01 (KM/SEC) *  IN  STATE 
TRANSITION  MATRIX 

VTOTPR 

A. 2b 

VELOCITY  (KM/SEC) * OUTPUT 

WASC 

SEE  INPUT  LISTING 

DEARTH 

17 

1 

RATE  OF  ROTATION  OF  THE  EARTH  (RAD/SEC)  *«*•>£ 
=•000072921150 

mIUD(7*2) 

BOO  INFORMATION 

mPERI 

SEE  INPUT  LISTING 

MX 

lbdA 

1 

INTERMEDIATE  VARIABLE 

60 


I 


I 


VARIABLE 

EOUA 

REF 

definition 

VAR  I 

rfXYZ 

169 

1 

INTERMEDIATE  VARIABLE 

XMU 

wY 

16BB 

1 

INTERMEDIATE  VARIARLE 

XN2  ( 

»  aBAR 

j 

89A  # 
97A 

1 

INTERMEDIATE  VARIABLE 

ANB 

AlG(brb) 

107 

1 

INTERMEDIATE  VARIABLES 

XWO 

AIMET 

INTERMEDIATE  VARIABLE 

AWUP 

XlNCHM 

INCREMENT  OF  ALTITUDE  IN  DNP  TABLE  =  10.0  KM 

AINCTM 

INCREMENT  OF  TEMPERATURE  IN  DNP  TABLE  =  -50.0 

DtG  K 

AW1 

XlPET 

INTERMEDIATE  VARIABLE 

AW1P 

XISCLH 

CONTANT  USED  TO  COMPUTE  POINTS  TO  BE  USED 

FOR  HIGH  ALTITUDE  EXPONENTIAL  EXTRAPOLATION 
=50.0 

XW2 

XW20 

A I  SUN 

177A 

1 

DIRECTION  COSINE  OF  SUN 

AW2P 

XLUEN 

INTERMEDIATE  VARIABLE 

XW2P 

XLST 

EXOSPHERIC  TEMPERATURE  BELOW  WHICH  THE  OENSITY 
WILL  NOT  BE  COMPUTED  =  500  DEG  K 

XWK 

AM(3»3) 

ASSIGNED  VARIABLES 

AWKP 

amean 

A. 57 

SEE  INPUT  LISTING 

AWYW 

AMINUN 

'19 

•31623*10  * 

XYSQ 

XMIPHT12) 

SEE  INPUT  LISTING 

AYZO 

XMN 

A. 19 

INTERMEDIATE  VARIABLE 

XYZS 

AMN1 

INTERMEDIATE  VARIABLE 

XYZS 

xmntru 

INTERMEDIATE  VARIABLE 

AYzT 

XMTR (3»  3) 

122- 

139 

1 

TRANSPOSE  OF  THE  XM  MATRIX  (SEE  XM.DEFINED 

IN  SUBROUTINE  TRANSF) 

61 


^IABlE 

EOUA 

REF 

DEFINITION 

J 

INTERMEDIATE 

VARIABLE 

2(3»3) 

14 

1 

NUTATION  MATRIX 

B 

3. o» Interval 
SYSTEM 

BETWEEN  1 

EPOCHS  OF  BASIC  REFERENCE 

U 

64»lu3» 

105 

1 

INTERMEDIATE 

VARIABLE 

(UNPERTURBED 

VALUE) 

UP 

b2» 103 » 
105 

1 

INTERMEDIATE 

VARIABLE 

(UNPERTURBED 

VALUE) 

1 

84 

1 

INTERMEDIATE  variable 
PERTURBED  VALUE) 

(UNPERTURBED 

OR 

IP 

b2»84 

1 

INTERMEDIATE  VARIABLE 
PERTURBED  VALUE) 

(UNPERTURBED 

OR 

2 

98A 

1 

INTERMEDIATE 

VARIABLE 

20 

98A* 103 

1 

INTERMEDIATE 

VARIABLE 

(UNPERTURBED 

VALUE) 

2P 

1U0A 

1 

INTERMEDIATE 

VARIABLE 

2P0 

lOuA  » 10b 

1 

INTERMEDIATE 

VARIABLE 

(UNPERTURBED 

VALUE) 

K 

A  *24 

INTERMEDIATE 

variable 

*P 

A. 2o 

INTERMEDIATE 

VARIABLE 

rwMU 

A. 26 

INTERMEDIATE 

VARIABLE 

bQKT 

INTERMEDIATE 

variable 

20(3) 

lb 

1 

SATELLITE  COORDINATES 

IN  CELESTIAL  SYSTEM 

2S<3) 

lb 

1 

STATION  COORDINATES 

2S0(3) 

126- 

133 

1 

STATION  COORDINATES  IN 

celestial  system 

2T(3) 

lo 

1 

SATELLITE  COORDINATES 

(X»  Y pZ)  IN  A  TOPOCEN- 

THIL  SYSTEM 
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VARIABLE 

EOUA 

REF 

definition 

AYZTP(O) 

17 

1 

SATELLITE  COORDINATES 
TRlC  SYSTEM 

(X»Y»Z)  IN  A 

TOPOCEN 

YBAR 

696 '97b 

1 

INTERMEDIATE  VARIABLE 

TWO 

83 

1 

INTERMEDIATE  VARIABLE 

(UNPERTURBED 

VALUE) 

YWOP  82' 103-106 

1 

INTERMEDIATE  VARIABLE 

(UNPERTURBED 

VALUE) 

YW1P 

62 '84 

1 

INTERMEDIATE  VARIABLE 
PERTURBED  VALUE) 

(UNPERTURBED 

OR 

YW2 

98b 

1 

INTERMEDIATE  VARIABLE 

YW20  986' 103' 104 

1 

INTERMEDIATE  VARIABLE 

(UNPERTURBED 

VALUE) 

YW2P 

1006 

1 

INTERMEDIATE  variable 

YM2P0  10GB' 105- 

106 

1 

INTERMEDIATE  VARIABLE 

(UNPERTURBED 

VALUE) 

YWKP 

A. 27 

INTERMEDIATE  variable 

ZDAT (3) 

SEE  INPUT  LISTING 

OR  HOUR 'MINUTE' SECOND 

OF  OBSERVATION 

ZETSUN 

177C 

1 

DIRECTION  COSINE  OF  SUN 

ZH 

ALTITUDE  OF  SATELLITE 

(KM) 

Z0NHAR15) 

156* 159 

1 

SEE  INPUT  LISTING 

ZR1 

156-159 

1 

INTERMEDIATE  VARIABLE 

ZK2 

158-159 

1 

INTERMEDIATE  VARIABLE 

ZR3 

158-159 

1 

INTERMEDIATE  VARIABLE 

ZR4 

158-159 

1 

INTERMEDIATE  VARIABLE 

ZRUXT 

134' 13o 

1 

INTERMEDIATE  VARIABLE 

ZSECU(3) 

149A- 

149C 

1 

ACCELERATIONS  X»V»Z  (KM/SEC2) 
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VI.  FLOWCHARTS 
A.  MAIN  PROGRAM  LOGIC 


MAIN? 


THIS  IS  THc  MAIN  PROGRAM  OF  THE  MAIN  LlNK.NOTE»THE  PROGRAM  IS 
LHAlhEo  AND  wRiTlEN  IN  FORTRAN  4. 


VARIABLE  EuUA  REF  DtF INiTION 

MERASE(I)  =U»FIKST  TIME  THROUGH  THE  PROGRAM  LOOP 

=1# MULTIPLE  PROBLEM 
OR 

=  1  *  ALL  DATA  HAS  BEEN  READ.NO  DISABLING 
DATA  ERRORS 

=  -1*  DATA  INPUT  ERROR. PROGRAM  WILL  GO  TO 
NEXT  PROBLEM 
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AFILT2 


THIS  IS  THt  MAIN  PROGRAM  OF  THE  FIRST  DEPENDENT  LINK • AFILT2  READS 
THE  UATA  (EXCEPT  FOR  THE  OBSERVATION  CARDS) .PRINTS  OUT  THE  INPUT  DATA. 
CHECtsS  FOR  Any  INPUT  DATA  ERROR  AND  COMPUTES  CERTAIN  INITIALIZATION. 

IF  THERE  ARE  ANY  DATA  ERRORS. THE  PROGRAM  WILL  PRINT  OUT  AN  APPROPRIATE 
ERROR  MESSAGE  ANu  WILL  CONTINUE  WITH  THE  NEXT  PROBLEM. 


VARIABLE  EQUA  REF 

LRASoU-2) 

LRASE ( 1*4 ) 

MERASEll) 


MERAbE(2-3) 


DEFINITION 

INTERMEDIATE  VARIABLES 

intermediate  variables 

=0. FIRST  TIME  THROUGH  THE  PROGRAM  LOOP 
=1 .NOT  THE  FIRST  TIME  THROUGH  THE  PROGRAM  LOOP 
OR 

StE  INPUT  LISTING 
OR 

=-l. INPUT  ERROR. PROGRAM  WILL  NOT  CONTINUE 
=1 .NO  OISARLInG  INPUT  ERROR .PROGRAM  WILL 
CONTINUE 

♦>  ' 

storage  cells 
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FIG.  6. 
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ERASE,  -  >7*1.0  4  (RASE |  *  (A. 20)  ♦  .00*2  *  (RASE,) 

CALL  JOATE  (I,  i.  MIS,,  ERASE,) 

CAU  JOATE  (I,  ),  MIS,,  ERASE,) 

ERASER  •  2*/)«S.O 

TREROS  •  (RASE,  «  (.»  ♦  .15  *  SIN  (ERASE^  *  (DJULN  -  ERASE,))) 
•  TNEXOS  *  SIN  (2.0  •  ERASE^  *  (DJULN  -  ERASE,)) 


ULL 

MIUR  (0) 


>sOITA  IN  MR  1 00" 


"CONFUTE  -  TWO  DAL  -  NODAL  MRIOO" 


(RASE,  -  ECCEN  •  ECCEN  ERASE%  ■  1 .0  -  ERASE, 


PNOOAL  .  ).0  •  ERASE ,/2.0  -  (SIN2  (OINCL)  *(*.0  •  .75  *  (RASE,) 
FNOOAL  -  MRIOO  *  (.  .0  -  1.5  *  TOMAN,  *  *  "« 

WRITE  (ON  IOUT)i  MRIOO,  ERASE.  TNEXOS,  PNQOAL  ' 


"READ  ALL  OISERVAT  ION 
CAROS  S  STORE  ON  *" 

— OsiD — 


TNORMN  -  OlDTIH 
DJULN  -  OLOOAV 


UFILT2 


THIS  IS  THE  ,v,AIN  PROGRAM  OK  THE  SECOND  DEPENDENT  LINK. THE  MAIN 

funciion  is  to  compute  ThE  Filtering  eouations*to  print  the  output# 

TO  RtJECT  ANY  OBSERVATION  CARDS  AND  TO  ITERATE. 


VARIABLE 

EUUA 

REF 

DEFINITION 

ol AG (6#6) 

52 

2 

1 

INTERMEDIATE  MATRIX  OR 

ThE  INVERSE  OF  THE  PHI  MATRIX 

lKASU ( 1 ) 

INTERMEDIATE  VARIABLE 

LRASc.U-20) 

INTERMEDIATE  variables  OR 

lRaSE ( 1*9) 

A. Id 

IMPROVED  STANDARD  OBSERVATION  ERRORS  COMPUTED 

AT  THt  END  OF  EACH  FILTERING  PHASE#OUTPUT 

lRASl( i-b) 

ly 

l 

CORRECTION  TO  POSITION  AND  VELOCITY  VECTOR# 
OUTPUT 

lRaSE ( 7-12) 

1 

STANDARD  DEVIATI0NS(SGUARE  ROOT  OF  THE  DIAGONAL 
ELEMENTS  OF  PMAT)  OF  THE  POSITION  AND  VELOCITY 
ERROR# OUTPUT 

tRASt ( 7 ) 

A  •  ltj 

INTERMEDIATE  VARIABLE 

lRaSE ( 19) 

A.lbA 

intermediate  variable 

cRASL  ( ) 

A.16A 

INTERMEDIATE  VARIABLE 

mPHTR ( 7  »  7) 
HPhTR (o»o) 

5b 

d 

1 

R  MATRIX#  OR  ITS  INVERSE  OR 

PHI  MATRIX  OR  ITS  INVERSE 

iCOUNT 

ITERATION  COUNTER 

MERASE(l) 

=-20# ALTITUDE  BELOW  120.0(KM)»  OR  EXOSPHERIC 
TEMPERATURE  OUTSIDE  OF  500aK  -  2400°K 

Range. PROgRAM  WILL  NOT  CONTINUE 
=i»COMPUTt  HIGHER  ORDER  TERMS  OF  THE  EARTH'S 
GRAVITATIONAL  POTENTIAL  FUNCTION  AND  THE 

DRAG  TERM  IN  THE  DIFFERENTIAL  EQUATION 
SUBROUTINE 

Pl2(b» J> 

lb 

1 

NUTAT lON-PRECtSSION  MATRIX  (COMPUTED  AT  DJULN) 

70 


H2  C  3»3) 
HHlTK(oro) 


13  1  PkECEbSION  MA'.MX  (COMPUTED  AT  DJULN) 

52  1  INTERMEDIATE  MATRIX 
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•  o 
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B.  SUBROUTINES 


bUBROUl INE  COATE  AND  SUBROUTINE  JDATE 

CDATE  CALCULATES  ThE  CALtNOAR  DATE  ( DAY » MONTH t YEAR )  GIVEN  THE 
MOolF ItD  JULIAN  DATE. 

uDATE  CALCULATES  THE  MODIFIED  JULIAN  DATE  GIVEN  THE  CALENDAR 
DATE ( DAY  t MON  TH  t Yt AR ) 


NOTE  1)  THE  MODIFIED  JULIAN  DATE  IS  THE  NUMBER  OF  INTEGRAL 
DAYb  SINCE  JANUARY  1»1950  (OhUT) 

k)  THE  YEAR  IS  DEFINED  TO  BE  THE  LAST  2  DIGITS  OF  19XX. 

3)  THE  VALID  CALENDAR  DATES  FOR  THIS  SUBROUTINE  ARE 

JANUARY  1 9 1950  (0hUT>  TO 

DECtMbEK  1  *  1999  C0KUT> 

4)  THE  VARIABLES  USED  IN  THE  2  SUBROUTINES  HAVE  NOT  BEEr 
DEFINED  IN  THE  LIST  OF  SYMBOLS. 


75 


ENTER  JDATE  \  (  ENTER  CDATE 

[MOAY,  MON,  MYR,  DAYJUL] )  V  [MDAY ,  MON,  MYR,  DAYJUL] 
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RETURN  TO 

CALLING 

PROGRAM 


SUBROUTINE  DENS IT 

THIS  SUBROUTINE  COMPUTES  THE  DENSITY  AS  A  FUNCTION  OF  ALTITUDE 
ANU  EXOSPHERIC  TEMPERATURE. 


VARIABLE 

EOUA 

REF 

DEFINITION 

ERASE (i) 

DISTANCE  OF  THE  SATELLITE  FROM  THE  CENTER  OF 
ThE  EARTH(KM) 

ERASE (20) 

b 

2 

EXOSPHERIC  TEMPERATURE  (DEG  K) 

TEMP 

177FW 

177C 

1 

INTERMEDIATE  VARIABLE 
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FIG.  12. 
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kJi 


< 


> 


U  -  1.0  +  (T-ETM)/X INCTM 
t  -  NTH- 2 


j  -  MAX0{ ( ZH  -  SHM)/X 1 NCHM,  ! } 

_ J 

l 

K  •  J  +  L  •  1 

DENl  -  INTERP  (DNP1K,  DNPU|  k,  DNP,+2  k) 

L  -  1.2,3 

1 

r 

U  -  (ZH  -  FLOAT  (ISH 
XLDEN  -  INTERP  (DENL 

DENTEM  -  (I0.0)XLDEN 

♦  j  *  1 NCHM) )/XI NCHM 

, ,  denl2,  denl3) 

INTERP  (ARG1 ,  ARG2,  ARG3)  - 

ARG2  +  275  *  j  ARG3  -  ARGl  +  U*  (ARG3  -  2.0  *  ARG2  +  ARG 1 )| 


FIG.  13- 
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SUBROUTINE  D1ACHK 

THE  MAIN  FUNCTION  OF  THIS  SUBROUTINE  IS  TO  COMPUTE  EQUATION  51 
REFEHENCt  1. 


VARIABLE 

EQUA 

REF 

definition 

0IAG(6'6) 

51 

1 

INTERMEDIATE 

MATRIX 

LRASO(l) 

INTERMEDIATE 

VARIABLE 

ERASE (1) 

INTERMEDIATE 

VARIABLE 

MERASE  ( 1 ) 

STORAGE  CELL 

PHITR (b»6) 

51 

1 

TRANSPOSE  OF 

THE  PHI (6*6)  MATRIX 

CALL  DIACHR  uuu, 


"*'1,1 

<  0.0 


"COVARIANCE 
MATRIX  IS 
NEGATIVE" 


CALL  SYXWT  (PMAT)J 


FIG.  14. 


1 


SUBROUTINE  EXROT 

THIS  SUBROUTINE  COMPUTES  THE  PRECESSION  MATRIX  (EQUATION  13*REF 
ERENCE  1.) 


82 


o 

CL 

X 


< 
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CL 

< 

H- 
to 

v_y 
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♦ 
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CM 
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H- 

CM 
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O 
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i 

— 
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O 

H 
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l 
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vO 

CM 

CO 

* 
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O 

00 

<A 

* 

CM 

o 

(A 

p— 
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CM 

CA 

• 

CA 
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• 

LA 

• 

CM 

p— 

00 

1 

CM 

1 

o 
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O 

CM  CA 
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— 

CA 
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• 
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-O 

• 

-O 

• 
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X 

X  X 

X  X 

X 

X  X 

LT\ 


O 


r\ 
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RETURN 


SUBROUTINE  HNATRX 

THIS  SUBROUTINE  COMPUTES  THE  DECLINATION  » RIGHT  ASCENSION  AND 
THE  h  MATRIX  (SECTION  Vll  BrKEFERENCE  1) 


VARIABLE 

EGUA 

REF 

DEFINITION 

CON (1-3) 

16*126* 

127 

1 

INTERMEDIATE 

VARIABLES 

CON (8) 

17*147 

1 

INTERMEDIATE 

VAR I ABLE *2  OR  Z 

ERASE (1-8) 

122- 

148 

1 

INTERMEDIATE 

VARIABLES 

MERASE(l) 

=l*STATION  CAN  SEE  THE  SATELLITE 
=-l*STATION  CAN  NOT  SEE  THE  SATELLITE 
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"COMPUTE :  Jf  Jf  # 


CALL  TRANSP(PNITR.  PNI ,3,3. 

ERASE,  -  -  COMj/(R  X  *  RT2) 

ERASE j  -  CON,* ERASE, 

ERASE,  -  COM,*ERASE, 

ERASE j  •  RX/RT2 

CALL  MULTIf  (ERASE,,.  PNITR,  ERASE.  3,  3.  I) 


1  _ —£ - 

M 

"COMPOTE :  ft-  K  , 

'  . 

ERASE,  •  -  CON2/RX2 
ERASE,  -  CON,/RX2 

CALL  MULTIP  ( ERASE j, 

H2 , 1  *  EMStl+2  !  1  ’ 

PNITR,  ERASE,  2.  2.  1) 
1.  2.  3 

RIP  .  100  . 


ERASE,  -  DV, 

>  1  -  1,  2.  3 

CALL  MULTIP 

(XT 20,  PNI,  ERASE,  3.  3.  ■) 

CALL  MULTIP 

(XY2S0,  PNI,  XY2S,  3,  3.  ■) 

CON,  -  XT 20 

-  XYZSO,  i  1  -  1,  2,  3 

RX2  -  CON,  *  COR,  ♦  CON,  *  CON,  '  ^COMPUTE  DECLINATION,  RIGHT  ASCENS-ON 
RX  -  (RX2)1'2  '*  *’ 


RT2  -  RX2  ♦  CONj  *  CONj 
OBSC,  -  AATAN  (CON,,  CON,) 
OBSC,  -  AATAN  (CONj,  RX) 

CALL  QUAO  (CON,.  CON,,  OBSC,) 


FIG.  16. 


MNRMHrOTM 


COMPUTE  I  ll ,  II .  H 


ERASEj  -  XYZTj/(RTC2  *  CONj) 

ERASE,  -  XYZT ,  *  ERASEj 

ERASEj  •  •  XYZTj  *  ERASEj 

ERASEj  •  C0N./RTC2 

CALL  MULTIP  (ERASEj,  XMTR,  ERASE,, 

3.  3,  1) 

H,  ,  •  ERASE, :  1  -  1.  2.  3 

'•COMPUTE  i  || ,  || ,  || 


ERASE,  •  XYZTj/CONj 
ERASEj  -  •  XYZT,/C0Nj 

CALL  MULTIP  (ERASEj,  XMTR,  ERASE,.  3,  2.  I) 
H2(|  •  ERASE|+2  i  I  -  I.  2.  3 


•R-  AR-  AR-  »Rt  »Rt  «Rt 

"COMPUTE  i  —  ,  —A,  — l-,  -rS  -4.  -r1 

•X  AY  12  *X  »Y  AZ 


COMPUTE:  4,  4.  4.  4 

•  X  AY  #2  «X  AY  (2 


ERASE,  •  XY2T j/CONr i  ERASEj 


OBSCj  ♦  OBSCj 


ERASEj  -  ((XY2TP3  *  XY2T ,  -  ERASE,  *  (XY2TP,  *  CONj  -  2R0XT  *  XYZT,))/CON5 


-  XV2T,  *  ERASEj J/RTC2 


ERASE,,  -  ( (XYZTPj  *  XYZTj  -  ERASE,  *  (XYZTPj  *  CONj  -  2R0XT  *  XYZTj ) )/CONj 
-  XYZTj  *  ERASE2)/RTC2 
ERASEj  -  (-  2ROXT  ♦  XYZTj  *  ERASEj )/RTC2 
CALL  MULTIP  (ERASEj,  XMTR,  ERASEj,  3,  3.  I) 


HS. 1*3  “  Ml , I 


I  -  I,  2.  3 


"COMPUTE:  4,  4,  4.  4 

_ _ _ i _ AX  AY  AZ  AX  AY  A  2 

ERASE,  •  OISCj  *  OBSCj 
ERASEj  .  -  (XYZTPj  ♦  ERASE,  *  XYZT,)/CONj 
ERASEj  -  (XY2TP,  -  ERASE,  *  XYZTj)/CONj 
CALL  MULTIP  (ERASEj,  XMTR,  ERASEj,  3.  2,  I) 

h6,  I  '  £MSEL3  |  , 

H6.U3-H3.I  r-'-2’3 


ERASE,  •  OBSCj/ATC2 
ERASE , +,  -  COM,  *  ERASE, 

Hj  ,  -  CON,tj/OBSCj  -  ERASE 

HA,  I  +3  *  MJ.  • 


"COMPUTE:  — 1,  — -4,  -X 

AX  AY  AZ  AX  AY  AZ 


j  I  -  I.  2,  3 


MTP  -  A 


COMPUTE:  -X.  “4. 

AX  AY _ AZ _ •  « _ *T _ A 1 _ *_ 


ERASE,  -  OBSCj  ♦  OBSCj  :  CONg  -  ZSECOj 
Hj  ,  -  (OBSCj  *  (CONj  -  ERASE,  *  Hj  , ) 
-  OBSC?  *  CON,  )/RTC2 
M7,UJ  *  2  0  *  (OVP,/OBSCj  -  ERASE, +  |) 


I  -  I,  2,  3 

j  .  II  -  I 


H.  ,  -  0,0 

|J:  A.  S.  6; 


FIG.  17. 
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SUBROUTINE  I NOBS 


THIS  SUBROUTINE  READS  ALL  THE  OBSERVATION  CARDS* CHECKS  FOR  ANY 
DATA  ERRORS  AND  STORES  ALL  ACCEPTED  CARDS  ON  A  DISK. 


VARIABLE 
ERASD(l) 
LRASEU-3) 
ERASE (b) 
ERASE (to) 
ERASE (7) 
ERASE (b) 
MERASE ( 1 ) 

MERASE(5> 


EQUA  REF  DEFINITION 

STORAGE  CELL 

INPUT  CELLS  FOR  AZIMUTH  OR  RIGHT  ASCENSION 

INPUT  CELL  FOR  RANGE  RATE 

INPUT  CELL  FOR  ELEVATION  RATE 

INPUT  CELL  FOR  AZIMUTH  RATE 

INPUT  CELL  FOR  RANGE  ACCELERATION 

INPUT  CELL  FOR  SATELLITE  NUMBER  OR 
=1* ALL  OBSERVATION  CARDS  MERE  NOT  IN  ERROR. 
=-l'ALL  OBSER VAT IONC ARDS  MERE  REJECTED  BE¬ 
CAUSE  OF  ERRORS. 

INPUT  CELL  FOR  ELEVATION  OR  DECLINATION 
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•  0 

1 s 1  uu  •  i  .71 


T  EM  I  .  ( 0  BSH |  ♦  ONUH  *  >0.0)  *  SIGMA 
OBM.  .  TCMI  *  FLR'.O 


"DEFINE  ELEVATION 
ON  DECLINATION" 


"DEFINE  RIGHT  ASCENSION" 


OBSMj  •  (ERASE,  ♦  ERASEj/60.0 
ERASE  ,/ 3600.0)  *  .26179938 

NTR  •  I 
NNTR  -  2 


FIG.  18. 
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OEFINE  AZIMUTH" 


FIG.  19. 
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bUBROUTlNE  INTEG 


THIS 

VARIABLE 
ERASE <C) 

tRASE(4-o) 

1COUNT 

MERASE(l) 

TEMP (6) 


lb  THE  INTEGRATION  SUBROUTINE 


EQUA 

REF 

definition 

155A 

1 

DISTANCE  OF  THE  SATELLITE  FROM  THE  CENTER  OF 
THE  EARTH  (KM) 

149A- 

149C 

1 

ACCELERATIONS  X » Y  r  Z  (KM/SEC2) 

ITERATION  COUNTER 

=-20  *  ALTITUDE  BELOW  120.0(KM)>  OR  EXOSPHERIC 
TEMPERATURE  OUTSIDE  OF  500°K  -  2400* K 
RANGE. PROGRAM  WILL  NOT  CONTINUE 
=1# COMPUTE  HIGHER  ORDER  TERMS  OF  THE  EARTH'S 
GRAVITATIONAL  POTENTIAL  FUNCTION  AND  THE 
DRAG  TERM  IN  THE  DIFFERENTIAL  EQUATION 
SUBROUTINE 

INTERMEDIATE  VARIABLES 


STM] 


CALL  INTCG 


r 

i 

i 


FIG.  20. 
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SU6R0U1 INE  1NVEKS 

THIS  SUBROUTINE 

VARIABLE  EUUA  REF 

tRASU(l) 

MERASE(l) 

HPHTR ( I » J) 
i=l#...»NNTP 
J=l» ...»NNTP 


INVERTS  A  MATRIX 

DEFINITION 

DETERMINANT  OF  HPHTR-  MATRIX  TO  BE  INVERTED 

RANK  OF  MATRIX  TO  BE  INVERTED 

MATRIX  TO  BE  INVERTED.  THE  INVERSE  OF  HPHTR 
IS  STORED  IN  HPHTR 
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FIG.  21. 
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SUBROUTINE  KEPLER 


1 HE  MAIN  FUNCTION  OF  Thlb  SUBROUTINE  lb  TO  COMPUTE  THE  PERIOD 
AND  THE  CONSTANT  K  (EQUATION  193#REFERENCE  1) 

VARIABLE  EOUA  REF  DEFINITION 

eRASE(I)  A. 31  COMPUTED  ECCENTRICITY  AT  THE  EPOCH  OF  THE 

INPUT  ORBITAL  ELEMENTS » OUTPUT 
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START  )  CALL  KEPLER  (KARG) 


«M 

CL  CM 
>  \ 
Q  — 


3  O. 

x  — 
o  -• 

\  UJ 
<j  UJ 

CM -I  * 


<■> 

CL  O 

CL< 

CM 

* 

CL 

—  I 

cn 

CM  CO 

5  * 

— J  to 

— 

CL 

CM 

>- 

CM 

UJ  w 

•  /• 

> 

> 

X  3 

(L 

UJ 

J-  3 

CM  0"\ 

O 

Q 

3E 

*  x 

— 

* 

•  Z 

> 

*  O 

CL  5 

—I 

l 

0)5 

Q 

+ 

+ 

N 

3iC>- 

UJ 

< 

CL  C 

ZJS 

QQ  CM 

o  -J 

■K  5 

+ 

CM  CM 

«— 

ZZ 

i  XX 

•  _J 

X 

CL 

Q. 

Si 

>- 

+ 

•— 

-< 

Cf\ 

CM  CM 

> 

> 

X  >- 

5*  * 

•* 

> 

O 

Q 

X 

CM 

CM 

^to 

CM 

O 

1  * 

3  3 

CL 

Q.  Q. 

CL 

r— 

+ 

* 

i  rr 

— 

—  — 

—  * 

+  X 

+ 

QbC 

— 1 

JJ 

-1 

h- 

CM  — 

, — 

o>:s 

UJ 

UJ  UJ 

UJ  t= 

O  Q 

CM  — 

CL 

> 

to  X 

< 

UJ  UJ 

O  CM 

•  v. 

> 

> 

Q 

> 

—  XX 

— 

Q  Q  w 

n  n 

i 

O  Q. 

Z  o’  Z. 

3c  </>  S 

X  >  X 


w  n 

1  3 

II  II  1 

II 

CM 

II 

II 

X 

CL  (L  CL 

CL 

CL  UJ 

cl.  3: 

^  — 

—  to 

-J  _l  _l 

-J 

-1  < 

Ul  UJ  UJ 

UJ 

UJ  CL 

>■  X 

<mu 

UJ 

UJ  UJ 

IUJ  I-  iC 
to  CL  O  o 


FIG.  22 
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SUBROUTINE  MOTION 

1HIS  SUBROUTINE  COMPUTES  THE  ECUATIONS  Oh  MOTION  (SECTION  Vlll  A 
VI 11  C .REFERENCE  1.  NOTE# THE  TERMS  INVOLVING  THE  TESSERAL  HARMONICS 


ARE  NOT 

USED  IN 

THE 

PROGRAM) 

VARIABLE 

EGUA 

REF 

definition 

ALTI 

SATELLITE  ALTITUDE  (KM) 

ERASE (A-tt) 

INTERMEDIATE  VARIABLES (EQUATIONS  OF  MOTION) 

ERASE(I) 

lb5A 

1 

ERASE ( 2 ) ♦ERASE ( 2 ) 

ERASE (l) 

1S5A 

1 

DISTANCE  OF  THE  SATELLITE  FROM 
ThE  EaRTH  (KM) 

THE  CENTER  OF 

ERASE (4— 6) 

149A- 

149C 

1 

ACCELERATIONS  X»Y»2  (KM/SEC2) 

MERASE(l) 

=-20  #  ALTITUDE  BELOW  120.0(KM)» 

OR  EXOSPHERIC 

TEMPERATURE  OUTSIDE  OF  500° K  -  2400* K 
RANGE. PROGRAM  WILL  NOT  CONTINUE 
=1 » COMPUTE  HIGHER  ORDER  TERMS  OF  THE  EARTH #S 
GRAVITATIONAL  POTENTIAL  FUNCTION  AND  THE 
DRAG  TERM  IN  THE  DIFFERENTIAL  EQUATION 
SUBROUTINE 
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ERASEj  -  DV,  *  ov, 
ERASEg  •  OVj  *  OVj 
ERASE,  -  C^ASCy  + 

ERASEj  -  (ERASE,)1 


♦  ov2  *  ov2 


♦  ERASE. 


Fig  -  -  39.375  *  ZR4  ♦  26.25  *  ZR2  -  1.875 
Flc  -  ZRI  *  (-  86.625  *  ZR4  ♦  78.75  *  2R2 
5  -  13.125) 

FI.  -  ZR2  *  (-  187.6875  *  ZR4  ♦  216.5625  * 
b  ZR2  -  59.0625)  ♦  2.1875 


FSMI  -  0.0 
FSM2  -  0.0 


CALL  MOTION 


ERASE UJ  ■  DV,  *  ERASE j 

I  -  I.  2.  3 


DVj/ERASEj 

ZRI  *  ZRI 
ZR2  *  ZRI 
-3.0  *  ZRI 

-7.5  *  ZR2  ♦  1.5 
-17.5  *  ZR3  ♦  7.5  *  ZRI 
ZR2  *  ZR2 


RER2  •  RER  *  RER2 

ZR2  -  RER2  *  ZONHAR, 

i  i  .  1 .  5 

FSMI  -  FSMI  ♦  ZRZ  *  r  1  |  ^ | , 
FSH2  -  FSH2  ♦  ZR2  *  F  1 , 

"POTENTIAL  FUNCTION 
IN  EQUATION  OF  MOTION" 


ZR2  -1.0+  FSMI 


ERASE,  -  ERASE,  * 

ZR2  :  I  •  4,  5.  6 

1 

* 

ERASEg  -  ERASEg  ♦ 

OMU  *  FSM2/ERASE, 

"COMPUTE  ALTITUDE' 
TN2PHI  •  ERASEg/ ERASEg 


WXYZ  -  (NX  *  WX  ♦  WY  *  WY  ♦  OVPj  *  OVPj)' 

ROPP  -  CORAG  *  F CONST  *  RHO  *  WXYZ 
ERASEg  -  ERASEg  +  WX  *  ROPP 

ERASEg  -  ERASEg  +  WY  *  ROPP 
ERASEg  -  ERASEg  ♦  OVPj  *  ROPP 


FIG.  23. 


SUBROUTINE  NUT ATX 


THIS  SUBKOUT XNE  COMPUTES  THE  NUTATION  MATRIX  (EQUATION  14* 
KEFEKENCE  1) 

VARIABLE  EUUA  REF  DEFINITION 

ERASE (1-10)  14  1  INTERMEDIATE  VARIABLES 
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i 


■4- 

CM 

l D 
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SUBROUTINE  ORBiNT 

THE  MAIN  FUNCTION  OF  THIS  SUBROUTINE  IS  TO  DETERMINE  THE  PRINT 
TIME  OF  THt  EPHEMERIS#  TO  INTERCONNECT  THE  FILTERING  AND  EPHEMERIS 
COMPUTATION#  AND  TO  tXECUTE  THL  SMOOTHING  OPTION. 


VARIABLE 

EQUA  REF 

DEFINITION 

UlAG(6»6) 

A  •  10 

INTERMEDIATE  MATRIX 

ERASE ( 1-12 ) 

INTERMEDIATE  VARIABLES 

ERASE (1-6) 

A. 11 

CORRECTION  IN  POSITION  AND  VELOCITY  DUE  TO 
SMOOTHING.  OR 

DIFFERENCES  BETWEEN  STORED  AND  INTEGRATED 
VALUES  IN  POSITION  AND  VELOCITY 

H(6»b) 

STORAGE  CELLS  USED  AS  A  MATRIX 

HPHTR (b»6) 

A.  10 » 

A. 11 

MATRIX  TO  BE  INVERTED  OR  ITS  INVERSE 

MERASE(l) 

=-20 .ALTITUDE  BELOW  120.0(KM>»  OR  EXOSPHERIC 
TEMPERATURE  OUTSIDE  OF  500®K  -  2400*K 
RANGE. PROGRAM  WILL  NOT  CONTINUE 
=1»C0MPUTE  HIGHER  ORDER  TERMS  OF  THE  EARTH'S 
GRAVITATIONAL  POTENTIAL  FUNCTION  AND  THE 
DRAG  TERM  IN  THE  DIFFERENTIAL  tOUATION 
SUBROUTINE 

bMATS(o»b) 

A. 11 

INVERSE  OF  THE  MATRIX 
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"DETERMINE  PROPER  PRINT  TIME" 


TIMREG  -  TIMREG  ♦  DPRINT 


rioolt  ft  U  XflMREG^X  >  J 

OJREG  -  DJREG  ♦  1 

|2°op  \86400 . 0^/  * 

TIMREG  -  TIMREG  -  86LOO.O 

KPIIPR 
•  0 


TNORMN  -  TSPI 
OJULN  -  DSP  I . 


JSP  •  JSP  *  I 
NOSPRI  •  NOSPRI  -  I 


IFILTR  -  I 


OJULN  •  ORYFLR 
TNORMN  .  TIMLR 


IFILTR 
•  0 


IFILTR  •  -I 


IFILTR  •  I 
TNORMN  •  T I  Ml  EG 
OJULN  •  DJREG 
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TNORMN  -  TIMFLA 
OJULN  -  OAYFLA 
CALL  DIACHK  (ICO) 
MTOTFLA  •  HTOTAL 


KFILPA  -  0 
I LMTO  -  I  LIMIT 
•LIMIT  -  0 


"SMOOTHING" 


ILF  IN  •  KIEG  ♦  I 


K  •  KIEG  NNTP  •  6 
OLDOAY  -  SMOAYk|[c 

OLDTIM  -  SMTIiy||(. 

WRITE  (ON  IOUT ) i 
OLDOAY,  OlOTIM 


"*T||,  y  "  *> 


SPPMT,  | .  ks  II  -  j 


""*ii.  j  •  "*Tn.  i  ♦  SPWTn.  k! 

II  -  1 . 6  J  -  1,....  6 

CALL  INVERS  NOlh  -  6 

SMATSUi  jj  -  HPHTA,,.  Jjt  II  -  I. 


caU  mult  ip” Jo?ag,  pmat,  SMATS,  6.  6,  6) 
HERASE,  -  I  CALL  INTEG 

ERASE,,  -  SOV,,  R  -  OVH,,i  II  -  I,...,  6 

WRITE  (ON  IOUT)! 

OJULN,  TNORMN,  ERASE  (1-6) 


CALL  MULTIP  (ERASE,  Dl AG,  ERASE,  6,  6,  I) 

OVH, ,  -  DVH, ,  ♦  ERASE  |  ut:  II  -  I . 6 

MPMTR|  | ,  j.  -  PHI,,  jj  i  II  -  1 . 6JJ  •  1 . 6 

CALL  INVERS 

PHI  1 1  ,ji  ■  MPHTR, ,  , ,  i  II  -  I . 6  JJ  -  I . 

CALL  FULTIP  (ERASE,  PHI,  ERASER,  6,  I) 


'I  I, KIEG 


S0VII,KIEG  *  tMStll 


DVII  "  S0VII,KIEG 
WRITE  (ON  IOUT)! 

ERASE  (1-6),  OV  (1-3),  DVP  (1-3) 


PMATT  II.  JJ  -  HI  I,  JJ 
II  —  I  ...o  J  ■  I ,  , . .6 


TNORMO  -  OLOTIM 
OJULO  •  OLDOAY 


I  UNTO 
•  0 


KANG  -  -I 
T  NO  AMI  •  OLDTIM 
OJULN  •  OLDOAY 


■  FILTR  -  0  I - HIS,. 


FIG.  27. 
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SUBROUTINE  P12SUB 


THIS  SUBROUTINE  COMPUTES  THE  TRANSFORMATION  MATRIX  BETWEEN 
TWO  BASIC  REFERENCE  SYSTEMS.  (SEE  EQUATION  PREFERENCE  1) 


VARIABLE 

EQUA 

REF 

definition 

P12(3#3) 

IS 

1 

TRANSFORMATION  MATRIX  OF  THE  RECTANGULAR  CO¬ 
ORDINATES  BETWEEN  TWO  SYSTEMS  OF  ARBITRARY 
DATES 

P2(3*3) 

13 

1 

PRECESSION  MATRIX 

(COMPUTED  AT  DJULN) 

105 


START 


CALL  PI 2  SUB 


FIG.  29. 


106 


SUBROUTINE  POSV/El 


1HIS  SUBROUTINE  CONVERTS  THE  CLASSICAL  ELEMENTS  TO  POSITION  AND 
VELOCITY  COORDINATES. (APPENDIX  6) 


VARIABLE  EOUA  REF 

ERASE ( 1 ) 

MERASE(l) 

H2 (3»3)  A. 56 

XM(2»2)  A • 54 » A. 5b 


DEFINITION 

INTERMEDIATE  VARIABLE 
ITERATION  COUNTER 
INTERMEDIATE  MATRIX 
INTERMEDIATE  MATRIX 
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START 


CALL  POSVEL  (XMEAN,  AXSEMI,  TSEC,  ECCEN,  WASC,  WPERI,  OINCL) 


SUdROUIlNE  PROUT 


THIS  SUBROUTINE  COMPUTES  THE  EPHEMERIS  AND  STORES  IT  ON  A  BCD 
ANU/OR  BINARY  TAPE. 


variable 

EUUA 

REF 

definition 

ALT  1 

A. 61 

NORMALIZED  SATELLITE  ALTITUDE 

LON(i-,i> 

16 

1 

INTERMEDIATE  VARIABLES 

iON(b) 

17 

1 

INTERMEDIATE  VARIABLE^ 

UJULO 

MODIFIED  JULIAN  DATE  OF  THE  PRINT  TIME 

LRASu( a-10) 

INTERMEDIATE  VARIABLES 

tRASuU) 

lu 

1 

GREENWICH  HOUR  ANGLE  OF  THE  VERNAL  EQUINOX 

LRASl>(J) 

11 

1 

INTERMEDIATE  VARIABLE 

LRASUUO) 

10 

1 

INTERMEDIATE  VARIABLE 

LRASE ( 1-20 ) 

INTERMEDIATE  VARIABLE 

tRASE(2) 

A. 24 

DISTANCE  OF  THE  SATELLITE  FROM  THE  CENTER 

OF  THE  EARTH  (KM) 

LRASE (b) 

A.  26 

VELOCITY  SQUARED C KM/SEC)2 

MERASE(l) 

STORAGE  CELL 

TNORMO 

TIME  OF  DAY  IN  SECONDS  CORRESPONDING  TO  DJULO 
(SEE  UJULO » ABOVE) 
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TNORHO  -  T IHSEC, 

UULO 

DJULO  -  DJUPRTj 

*  < 

s^OJUPUT 

II  .  LUf 

V' 


FIG.  31. 


FIG.  32. 


1 1  1 


POSITION-VELOCITY  COMPUTATION' 


ERA$E7  -  OV,  *  DV,  ♦  ov2  *  dv2 
ERASEg  -  DVj  *  OV3 


TN2PHI 


ERASEg/ ERASEg 


ERASEj  -  <ERASE7  +  ERASEg)172  RAOPRT  .  ERASE 


ERASEg  -  DVP2  ♦  OVP^  ♦  OVP^ 


ALT  10  «  ERASE. 


REARTH* 


VTOTPR  .  ( ERASEg ) 
1/2 


1/2 


ALTI  -  ALTIO/REARTH 


RADIUS 

VELOCITY  (TOTAL) 
ALTITUDE 


ERASE 


7 


(erase7) 


1/2 


GEOCEN  -  AATAN  (DVj,  ERASER 


SOGA  - 

SOGBO 
SOGB  • 


dv3/erase7 


•  SOGA 
SOGA 


I TCT  -  0 


GEOCENTRIC  LATITUDE 


"ITERATE  FOR  GEODETIC 
LATITUDE" 


GEOCEN  -  AATAN  (OVj,  0.0) 
GEODET  .  GEOCEN 


,nr,  /I.O  ♦  SOGB  *  SOGB 

vl/2 

-0GC  "^1.0  ♦  FEARTH  *  SOGB 
'OCD  /ALT!  ♦  S09?  1 

*  s6gb ) 

|  *  SOGA 

„OCO  l ALT j  ♦  FEARTH  *  SOGC  J 

1  TCT  -  ITCT  ♦  1 

sog^yom 


I 


<1.0*  10 

7 - 


GEODET 


TAN*1  (SOGB) 


GEODETIC  LATITUDE 


OLAMS  -  AATAN  (DV2>  DV , ) 
CALL  QUAD  (DV,  ,  DVj,  OLAMS) 


6  ■  V  1  ^ 

"GREENWICH  HOUR  ANGLE"  ^ 

' 

ERASDj  -  -.7665486I*I0'4*  (SIN  (.21140824  -  .924221651  *  I0"3  *  DJULO) 

-.567232*I0'5*  (SIN  (9.7766782  +  .03440557  *  DJULO) 

ERASD|0  -  1.746647719  +  .0I720279I45097297*DJUL0 

♦  .5064  *  I0'14  *  DJULO  *  DJULO  +  .000072921158  *  TNORMO 

ERASD  -  ERASD, Q  ♦  ERASDj 

"SATELLITE  LONGITUDE' 


"GREENWICH  MEAN  SIDEREAL  TIME 
TO  -  HOUR.  MINUTE,  SECONDS"  , 


OLAMS  -  ERASD^  -  OLAMS 
10  " 

^OLAMSP  -  2t  *  FLOAT  (MERASE ,  g) 


OLAMO  - 

GEODET  -  GEODET  'FLRAD 
GEOCEN  .  GEOCEN/FLRAD 


FLRAD 


-) 


MERASE 

erasd6 

ERASD, 


10 


ERASD,g/2t 
MERASE 


ERASD, 

MSTHR 

ERASD 


10 

* 

ERASD 


ERASDg  *  2* 


10 

MSTHR  -  ERASD 


10  •  ERASD$ 


2  » 


*  24.0 


ERASD|q 

MSTMIN 

ERASDg 

ERASD, 

STSECO 


10 
MSTHR 

■  ERASD, j  -  ERASDg 
ERASD, ,,  *  60.0 

FLOAT  (  MSTMIN  )  /  60.0 

■  ERASD..  -  ERASD, 


ERASD 


10 


3600.0 


FIG.  33. 
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EAASEg  -  (OV,  *  OVP,  ♦  OVj  *  OVPj  ♦  OVj  *  OVPjJ/ERASEj 

ERASEj  -  {ERASE,:  -  EAASEg  *  EAASEg)l/2 
EAASEg  •  ERASE  *  ERASEj/OHU 
EAASEg  ■  1.0  -  EAASEg  *  ERASE 7 

ERASE |0  -  ERASEg  *  ERASEg  ERASE,,  -  (ERASEg)2  *  OhU 
ERASE, 2  .  ERASEg  *  ERASE,  ♦  ERASE, „  *  ERASE, Q 

ECCENT  -  (ERASE, 2),/2  AXIHAJ  •  ERASE,, /( 1.0  •  ERASE, 2) 
ERASE, 2  -  (1.0  -  ERASE, 2),/2 

Aj  .  0V2  *  DVPj  -  OVj  *  OVPj  Aj  -  DVj  *  OVP,  -  DV,  *  OVPj 

AL  -  OV,  *  OVP,  -  OV,  *  OVP,  ERASE,,  -  (A}2  ♦  Aj2  ♦  Ag2)17 


Ag  -  OV,  *  OVPj  -  0V2  *  OVP, 
*1+3  "  *|/£RAStU  1 


A,  -  OV , . ,/ERASE 2 


"RIGHT  ASCENSION  OF  ASCENDING  NODE" 


ERASE, g  -  A,  ERASE, 7  -  -Ag 
ASCNOD  •  AATAN  (ERASE, g,  ERASE , 7 ) 
CALI  QUAD  (ERASE, 7,  ERASE,g,  ASCNOD) 
ASCNOD  -  ASCNOD/FLRAO 


i  -  I.  2,  3 


ERASE, g  •  (1.0  -  A7  *  A7)  inns 

OINCLI  -  AATAN  (ERASE, g.  ERASE, j) 
CALL  QUAD  (ERASE, j.  ERASE, g.  OINCLI) 
OINCLI  •  OINCLI/FLRAD 


ERASE, j  -  A j 


ERASE, g  -  Aj  *  ERASE,  7 


♦  *3  *  £**«|6 


ERASE, j  -  *;  *  (Aj  *  ERASE, 7  -  Aj  *  ERASE, g)  ♦  *g  *  (*5  *  CRASE ,6  -  Ag  *  ERASE ,7 ) 

ERASE, g  -  TAN*1  (ERASE , j/ERASE ,g) 

CALL  QUAD  (ERASE, g.  ERASE, j,  ERASE, g) 

ERASE, g  -  -ERASE, 0/ECCENT  ERASE, ?  -  -ERASEg/ECCENT 
ERASE, g  •  AATAN  ( ERASE, g,  ERASE , 7 ) 

CALL  QUAD  (ERASE , ERASE ,g,  ERASE, g) 

PERIGE  -  ERASE, g  ♦  ERASE, g  HERASE ,Q  •  PERIGE/2* 


PER  I  GE/  2  v 


PERIGE  *  (PERIGE  -  2*  *  FLOAT  (HERASE ln) )/FLRAO 


,  2 

♦  ECCENT 


ECANOM  -  AATAN  ( ERASE, g.  ERASE, j) 
CALL  QUAD  (ERASE, ERASE, g.  ECANOH) 
0  HE  ANA  .  ECANOH  -  ECCENT  »  ERASE,. 


OHEANA  -  OHEAN  ♦  2» 


OHEANA  .  OHEANA/FLRAD 


LFO  -  5 


LPAGE  •  LPAGE  ♦  1 
LINOUT  -  3 

*  >  < 

■^UNOUT  N- 

WRITE  PROPER 
HEADING 

<5^ 

M3 

LFO  -  7 


1/  KEVTAP 
\  <0 


LOO  AAR  -  0 


rttprt 
•  s 


"BCD  TAPE-EPHEHER I  $  OUTPUT" 

LOO  KAN  .  I  LINOUT  -  LINOUT  ♦  4 

WRITE  (ON  HSCD) :  OAYJL.  KHOOUT  KOAOUT ,  RYROUT ,  KHROUT, 
KMIOUT ,  SECOUT,  TNORM0,  IREV,  OV ( I  -  3 ) .  DVP(l-3), 

ALTIO,  RAOPRT,  VTOTPR,  GEOCEN,  GEODET,  OLAHO. 

HSTHR.  MSTHIN.  STSECO 


tft POLO  •  OTP 
Ij  .  I 
OTP  •  IOO 
OLAHS  -  ERASE  * 

I  •  0 


1.1*1 

OS  -  I 

ERASE  |  •  SLOO, 

RALOO  •  OLAMS  -  ERASD, 
CAU  HMATRX 


RTASC  •  OISCj/FLAAO  OECLIO  •  OOSC,/FLRAO 
CALL  RFCOR 

ORSM,/FLRAD  A2IHUT  .  OISH,/FLAAD 
OOSH,  HAORA T  -  01 SO^ 


AERATE  -  ORSCj/FLRAO 


ELRATE  -  OISCj/FLRAO 


ERASE^  -  COO,  •  COO,  *  COO,  *  COO, 

RARATE  >  (COO,  •  COO,  •  COO,  *  C00j)/(ERASC^  •  FLAAO) 

{CON,  •  ERASE.  -  COO.  *  (COO,  •  CON.  «  COO,  •  COO.) 

— 3 - 5 - 1 - —J - 1— - 1 - 1— 

RTCE  «  FLRAO  >'  (tRASE^) 


IJ  -  li  - 


LOO  KAN  .  I  LINOUT  •  LINOUT  *  2 

WRITE  (ON  MCO):  DAYJL,  WOOUT,  ADAOUT, 

KYROOT,  RHROUT,  AN  I  OUT,  SECOUT,  TOORMO,  IREV 


STATION  DATA  •  ICO  TAPE  -  OUTPUT 


LIOOUT  •  LINOUT  *  2 

WRITE  (ON  MCD):  NUHSTA^j.  ELRATE,  AERATE,  RARATE. 
DCRATE,  ELEVAT,  A2IMUT,  RANGES,  RANRAT,  RTASC.  OECLIO 


IJ  <  10 


Ij  -  10 

WORDS  .III 


NUMER,.  ITYPE 

HERASE,  •  NWOROS  *  1 

WRITE  (ON  MIN)  i 

HERASE,.  IJ, 

AHCOATIt  1  •  1.. 

..  MAUDS 

_ : _  .  _ i 

Ull  *  (IJ-I)  ♦  I 
ERj  +  ,  •  NUHSTAn 


"WRITE  EPHENERIS  ON 
-1  UNARY  TAPE" 


HERASE,  .  NWOROS  *  31  NUMER,  -  I  TYPE 
WRITE  (ON  MIN)i  HERASE , ,  DAYJL,  RMOUT,  KDAOUT, 

KYROUT,  RHROUT,  KMIOUT,  SECOUT,  TNOROS,  D¥(l-3).  DVP(l-3). 
ALT  10,  RADPRT ,  VTOTPR,  GEOCEN.  GEODET,  OLAHD.  MSTHR, 
HSTHIN,  STSECO,  AXIHAJ.  ECCENT.  OINCLI,  ASCOOD,  PERIGE, 
OMEANA,  IREV,  IJ.  (AOGOAT,!  j  .  I . NWOROS) 


I  TYPERS.  > 


ANGDATj  +  j  -  DCRATE 
AOGOATj  ^  -  A2IMUT 

ANGDATj  ,  9  .  RANRAT 
AOGOATj  ,  ,, 


ANGDATj  +  6  .  ELEYAT 
ANGDATj  4  j  -  RANGES 
AOGOAT,  .  -  RTASC 


IJ  •  IJ  • 


Ij  •  II  «  I 


.NEC  -  IREC  - 


OTP  -  NT  POLO 


LUP  •  I  SAVE 


II  -  II  «  I 


iSnacfil 

tlr'.fi 


SUBROUTINE  RFCOR 

THIS  SUBROUTINE  COMPUTES  THE  REFRACTION  CORRECTIONS. (SEE 


SECTION  IV 

D» REFERENCE 

1) 

VARIABLE 

EftUA  REF 

DEFINITION 

ERASE (1-20) 

INTERMEDIATE 

VARIABLES 

TEMPE 

INTERMEDIATE 

VARIABLE 

115 


CALL  RFCOR 


FIG.  36. 
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H7 


bUbROUTlNE  SETCON 


THIS  SUBROUTINE  COMPUTES  PROGRAM  CONSTANTS  AND  CERTAIN  INITIAL- 


IZATiON 

FROM  THE 

INPUT 

DATA. 

VARIABLE 

EOUA 

REF 

DEFINITION 

MERAbE ( 1 ) 

SEE  INPUT  LISTING 

ERASE (1-5) 

5-7 

1 

INTERMEDIATE  VARIABLES 

118 


CALL  SETCON 


FIG.  38. 
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SUBROUTINE  TRAMAT 

THIS  SUBROUTINE  COMPUTES  THE  STATE  TRANSITION  MATRIX  (SECTION 
VI  REFERENCE  1) 


VARIABLE 

EuUA 

REF 

definition 

PHITR(b'6) 

INTERMEDIATE 

MATRIX 

TEMP 

INTERMEDIATE 

variable 

1EMPE 

9b 

1 

INTERMEDIATE 

variable 

CALL  TRANAT 


0* 


SOS  -  XWI  *  YWIP/OMU 
SOC  .  OMU  *  SOS  *  SOS 
SOA  -  I  .0  -  SOS  *  YWIP 
SOB  •  SOS  *  XWIP 

ECCJ  -  SOA2  ♦  SOB2 

ECC  -  CECC2) 1/2 
SO  .  -  (1.0  -  ECC2 ) 

SOKB  -  (SOK)l/2 
SOKV  -  1.0/ SOS 

S0KC  * 


_ X _ 

SAPA  .  -1.0/ SOC 
SIPA  .  0.0 
ECA  -  HTOTAS/SOKC 
ECAO  -  ECA 


* 


1  XWI  .  (OV,2  ♦  DV  2 

♦0V,2)'/2 

XIG, 

1  “ 

DV,/XWI 

*ig4  4  - 

me, 

X  1C, 

1  ■ 

OVj/XWI 

XIGj  4  - 

XIC2. 

*,c 

• 

1 

DV,/XWI 

x,gG .4  ■ 

HIGj 

SAI  -  DYj 

*  OVP,  - 

OV,  *  OVP, 

SA2  -  OV 

*  OVP,  - 

OV,  *  OVPj 

SA3  .  OV. 

*  OVPj  - 

dv2  *  OVP, 

SO  • 

(SAI 

2  ♦  SA22 

*  SA32)1'2 

X.G, 

3  ‘ 

SAI /SO 

x  1  gA  ,  6  ‘ 

XIGI  .3 

XIG, 

3  * 

SA2/S0 

X,GS.6  ' 

XIG2.3 

X 1 G  3 

3  * 

SA3/S0 

x,G6.S  * 

X,G3.3 

XIG, 

2  * 

X,C2.3  * 

XIC3. 1  '  X,G3.3  * 

XIC2, 1 

XIG, 

2  ‘ 

XIC3.3  * 

XIGI,I  '  X,GI.3  * 

X,63.l 

xig4 

5  ' 

XICI.2 

X ICS . 5  ‘ 

X,G2.2 

*IC3J 

2  ‘ 

x,cl ,3  * 

XIG2.I  *  XIG2.3  * 

XIGI.I 

x.g6 

S  * 

X,C3.2 

XWIP 

•  OVP,  *  XIG, 

.1  *  0VP2  * 

X,G2. 

♦  OVP,  *  XIG,  , 

YWIP 

-  OVP,  *  XIG, 

.2  ♦  0V'2  * 

XIG2.i 

.  OVP,  *  XIG,  2 

XWO  -  XWI 

YWO 

-  0.0 

|  XWOP 

•  XWIP  YV*P 

-  YWIP 

_ X _ 

STONE  •  COS  (ECAD) 

STOU  •  1.0  -  STONE  *  STONE 

STOUSQD  -  (STOU) 1/2 
XBAR  .  (STONE  -  ECC)/SOK 
YBAA  .  STOUSQ/SOKB 


FIG.  39. 
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RETURN 


) 


FIG.  40. 
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SUBROUTINE  TRANSF 

THE  MAIN  FUNCTION  OF  THIS  SUBROUTINE  IS  TO  COMPUTE  THE  ELEVATION* 
AZIMUTH*RANGE*KANGE  kate»elevation  rate»azimuth  rate  ano  RANGE  ACCE- 


LERATION. 

(SECTION 

VII 

B*KEFERENCE  1) 

VARIABLE 

EQUA 

REF 

DEFINITION 

CON (1-4) 

16 

1 

INTERMEDIATE  VARIABLES 

CON (6) 

17 

1 

INTERMEDIATE  VARIABLE*Z 

MERASE(l) 

=1*  STATION  CAN  SEE  THE 
=-l*STATION  CAN  NOT  SEE 

satellite 

THE  SATELLITE 

XM ( 3  *  3 ) 

16* 

17 

1 

intermediate  variables 
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CALL  TftANSF 


•  <• 
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SUBROUTINE  AUDIT!  OR  SUBROUTINE  SUBTRA 

THIS  SUBROUTINE  ADDS  OR  SUBTRACTS  2  MATRICES* 


SUBROUTINE  MULTIP 

THIS  SUBROUTINE  MULTIPIES  2  MATRICES* (NOTEtIT  IS  ASSUMED  THAT 
THE  FORTRAN  DIMENSION  OF  EACH  MATRIX  IS  EQUAL) 


SUBROUTINE  TRANSP 

THIS  SUBROUTINE  OBTAINS  THE  TRANSPOSE  OF  A  MATRIX* 


SUBROUTINE  A AT AN 

THIS  SUBROUTINE  COMPUTES  THE  ARCTANGENT. (GIVEN  THE  SINE  AND 
COSINE) 


SUBROUTINE  QUAD 

THIS  SUBROUTINE  OBTAINS  THE  PROPER  QUADRANT  OF  AN  ANGLE  AFTER 
IT  HAS  BEEN  OBTAINED  BY  THE  ARCTANGENT  SUBROUTINE. 
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EX 

M 

) 

IBBBflysl 

-  1,...,  NCB } 

CALL  ADDITI  (A.B.C.NRB.NCB)  OR 
CALL  SUBTRA 


(A.B.C.NRB.NCB) 

(A.B.C.NRB.NCB) 


CALL  MULTIP  (A.B.C.NRB.NCB.NCC) 


CALL  TRANSP  (A.B.NRB.NCB) 
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SUBROUTINE  BCOFL  OK  SUBROUTINE  BCDFX 

SUBROUTINE  BCOFL  CONVERTS  A  BCD  WORD  TO  A  REAL (FLOATING  POINT) 
NUMBER •  SUBROUTINE  BCOFX  CONVERTS  A  BCD  WORD  TO  AN  INTEGER (FIXED 
POINT)  NUMBER.  THE  FORTRAN  CALLING  SEQUENCE  IS 

CALL  BCUFL ( ARG1 » N » ARG2 » J ) 

OR 

CALL  BCDFX ( ARG1 » N  > ARG2 * J ) 

ARG1-  WORD  TO  BE  CONVERTED 

N-  NUMBER  OF  CHARACTER  OF  ARG1  TO  BE  CONSIDERED  IN  THE 
CONVERSION 

AKG2-  THE  CONVERTED  NUMBER  (REAL  OR  INTEGER) t ALWAYS  POSITIVE. 
J-  =- 1» ERROR  WORD  CAN  NOT  BE  CONVERTED 
=0'  WORD  HAS  BEEN  CONVERTED 

NOTE.  THE  WORD  MUST  CONSISTS  OF  NUMERICAL  BCD  CHARACTERS .THERE  IS  NO 
FLOWCHART  AND  VARIABLES  OF  THE  SUBROUTINE  HAVE  NOT  BEEN  DEFINED. 
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SUBROUTINE  BCDSTO 


THIS  SUBROUTINE  IS  USED  TO  STORE  BCD  INFORMATION .THERE  IS  NO 
FLOWCHART. 


SUBROUTINE  SYMMET 

THIS  SUBROUTINE  SYMMETRIZES  THE  COVARIANCE  MATRIX  (EQUATION  51 
REFERENCE  1)  AND  THE  MATRIX  K*H*P  (EQUATION  52  REFERENCE  1.) 

THERE  IS  NO  FLOWCHAkt. 


SUBROUTINE  OVPUN  AimD  SUBROUTINE  NOVPU 

SUBROUTINE  OVPUN  CONVERTS  A  BCD  CHARACTER  TO  A  REAL (FLOATING 
POINT)  NUMBER.  SUBROUTINE  NOVPU  CONVERTS  A  BCD  CHARACTER  TO  AN 
INTEGER (FIXED  POINT)  NUMBER.  THE  FORTRAN  CALLING  SEQUENCE  IS 

CALL  OVPUN ( ARG1 » ARG2 » J ) 

OR 

CALL  NOVPU (ARG1 # ARG2» J) 

ARG1-  CHARACTER  TO  BE  INVERTED 

ARG2-  THE  CONVERTED  NUMBER  (REAL  OR  INTEGER) 

J-  =-l»ERFcOR. CHARACTER  CAN  NOT  BE  CONVERTED 
=0»  ARG2  IS  POSITIVE 
=1'  ARG2  IS  NEGATIVE 

THE  NET  RESULT  OF  THESE  SUBROUTINES  IS  TO  CONVERT  AN  OVERPUNCH (ON 
OBSERVATION  CARO)  TO  A  NUMBER .THERE  IS  NO  FLOWCHART  AND  VARIABLES 
OF  THE  SUBROUTINES  HAVE  NOT  BEEN  DEFINED. 
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APPENDIX 


A.  TIMING  ERROR 


A  timing  error  introduces  errors  in  satellite  observation,  which, 
obviously,  will  appear  along  the  direction  of  the  satellite  velocity  vector. 
Thus,  they  are  distinctly  different  from  normal  observation  errors. 


If  at  the  time  of  an  observation,  the  components  of  the  satellite  ve¬ 
locity  vector  are  x,  y,  z  and  the  components  of  the  acceleration  vector 
are  x,  y,  z,  the  error  in  position  due  to  a  standard  timing  error  a  will  be. 


Axt  =  Ixal  (A.  la) 

Ayt  =  |ya|  (A.  lb) 

Azt  =  |zol  (A.  lc) 

and  the  error  in  velocity  will  be 

Axt  =  |xal  (A.  2a) 

Ayt  =  |ya|  (A.  2b) 

Azt  =  |za|  (A.  2c) 


which  will  give  a  covariance  matrix  of  the  position  and  velocity  vector 
due  to  the  timing  error: 


Thus,  the  total  covariance  matrix  at  an  observation  time  including  the 
timing  error  will  be  (see  Ref.  1,  Section  V-D) 
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(A.  4) 


R  =  H  <tk+1>  *  (tk+1.  tk)  P'  <tk)  *T(tk+1,  tk)  HT  (tk+1) 

+  H  <W  pt  (tk+i>  »T  <W  +  «  <W 

which  can  be  rearranged 

R  =  H  (tk+l* [*  (tk+l*  V  P'  *T  (tk+l’  V 

+  Pt  'Vl']  H?  W  +  Q  <W  <A-5> 

This  covariance  matrix  is  expressed  in  terms  of  observations  and  re- 
places  the  matrix  inside  the  square  brackets  in  Eq  50,  Ref.  1. 

It  must  be  pointed  out  that  the  timing  error  cr  can  be  either  for  an 
individual  observation  or  a  particular  station  or  the  whole  system.  In 
a  general  case,  when  the  timing  error  for  a  particular  station  will  not 
be  known,  the  system  standard  timing  error  can  be  used  which,  gen¬ 
erally,  will  be  more  easily  estimated.  The  present  program  utilizes 
a  system  timing  error  which  is  an  input  value. 

B.  SMOOTHING  EQUATIONS 

Smoothing  can  be  considered  a  special  case  of  filtering  in  the  sense 
that  the  observations  are  replaced  by  the  state  variables  which  in  our 
case  are  the  six  orbital  elements  in  the  form  of  position  and  velocity 
vectors.  Returning  to  Eqs  (49),  (50),  (51)  and  (52)  of  Ref.  1: 


*  <w 

=  K 

<w 

[y'  <W 

■  y  <w] 

(A.  6) 

K  <w 

=  P 

^k+l* 

»T  <W  [ 

H<WP 

<W  «T  (tk+i> 

+ 

« <w] 

(A.  7) 

p  <w 

=  » 

(tk+r 

V  p'  <V. 

)  *T  (*15+1- 

*k> 

(A.  8) 

p,<w 

=  P 

(tk+l* 

- K  <w H  <w p  <w 

(A.  9) 

We  can  consider  that  instead  of  a  series  of  observations  obtained  at 
discrete  times,  we  have  a  series  of  estimates  of  the  state  variables 
and  their  covariance  matrices.  These  estimates  of  the  state  variables 
and  the  corresponding  covariance  matrices  are  obtained  in  the  process 
of  regular  filtering  and  stored  in  the  computer. 
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Assuming  that  at  time  t  we  have  an  estimate  x(t  It  )  based  on  ob- 

servations  obtained  up  to  and  including  time  tn,  we  can  update  a  previous 

estimate  at  time  t^,  x(t^|  based  on  data  which  includes  time  tR, 

(tk<  tn).  The  updated  estimate  will  be  designated  ^(t^l^).  The  process 

will  now  consist  of  computing  the  weighting  matrices  from  the  covari¬ 
ance  matrices  of  two  sets  of  estimates  at  each  data  point,  and  cor¬ 
recting  one  set  of  the  estimates. 

Returning  to  Eqs  (A.  6),  (A.  7),  (A.  8)  and  (A.  9)  we  can  make  the 
following  deductions.  Since  the  observations  are  now  replaced  by  es¬ 
timates  of  the  state  variables  or  the  position  and  velocity  vector,  the 
matrix  H  and  its  transpose  will  be  a  unit  matrix  because  the  compo¬ 
nents  of  the  position  and  velocity  vector  are  independent  variables. 

The  Q  matrix  will  now  be  the  covariance  matrix  P(tjJtk)  of  the  esti¬ 
mate  xit^lt^).  Then  we  can  obtain  a  weighting  matrix 

K  <tk)  «  P  (tk|tn)  [P  (tk|tn)  +  P  (tkl  tk>]  -1  (A.  10) 

The  updated  estimate  at  tfc  based  on  data  up  to  and  including  tn  will 
now  be 

5  «k'  y  ■ K  <v  [ 5  (tk'  v  -  *  <v  v  *  v]  (A- u> 

and  the  recursion  equation  for  the  updating  of  the  covariance  matrix 
will  be 

p'  <tk'V  =[l-K  <v]  •  <V  V  p  (‘n'V  *T  <V  V  <A- 12> 

In  the  present  application  of  the  smoothing  technique,  we  are  interested 
in  obtaining  the  best  estimate  of  the  orbit  at  time  tn.  Therefore,  after 

each  smoothing  operation  we  do  not  update  the  orbit  at  time  t^,  but  trans¬ 
fer  the  corrections  to  time  tn.  Then,  with  the  improved  orbit  we  pro¬ 
ceed  to  time  t^.^  and  repeat  the  smoothing  until  all  points  are  processed. 

Since  this  type  of  smoothing  process  may  involve  long  time  arcs  and 
cause  numerical  difficulties,  only  the  diagonal  elements  of  the  covari¬ 
ance  matrices  are  stored  during  the  filtering  process.  Similarly,  only 
the  diagonal  elements  of  the  P  (tn  I  tn)  matrix  are  used  in  the  transforma¬ 
tion.  In  addition  to  avoiding  numerical  difficulties,  this  approximation 
saves  considerable  storage  space.  However,  it  must  be  pointed 
out  that  the  P  (t^l  t  )  matrix  is  not  a  diagonal  matrix. 
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C.  SYSTEM  STANDARD  OBSERVATION  ERROR 


One  of  the  problems  in  a  sequential  orbit  estimation  technique  is 
the  estimation  of  standard  deviations  of  the  observations.  This  is 
somewhat  complicated,  as  discussed  previously  in  this  report,  because 
the  standard  deviations  of  the  system  measurement  error  include  un¬ 
known  bias  errors  as  well  as  random  errors.  A  solution  to  this 
problem  can  be  obtained  by  a  method  of  successive  approximations  in 
a  multiple  filtering  process. 

An  approximation  of  the  expected  deviations  of  the  measurements 
from  the  estimated  orbit  at  a  time  t^  can  be  obtained  from  the  diagonal 

elements  of  the  covariance  matrix. 

R  *  H  <tk)  P  (tk)  HT  <tk)  +  Q  <tk)  (A.  13) 

which  is  expressed  in  terms  of  the  observations. 


If  AS^  is  the  deviation  of  the  kth  measurement  from  the  estimated 
orbit,  then  we  can  write 


1.25 

n 


(A.  14) 


where  n  is  the  number  of  observations  and  r^  the  diagonal  element 

corresponding  to  the  particular  type  of  observation  in  matrix  R.  Both 

AS.^  and  ^r.^  include  the  deviations  in  the  estimated  orbit  and  the 

measurement.  In  the  case  of  normal  distribution  of  errors  and  if  the 
estimate  of  the  standard  deviation  is  correct,  then  should  approach 

1  for  a  sufficiently  large  number  of  observations.  If  Cj  4  1*  the  dis¬ 
crepancy  will  be  due  to  the  incorrect  estimates  of  the  standard  devia¬ 
tions  in  the  observations.  Therefore,  a  better  approximation  can  be 
obtained  assuming  that  the  error  in  the  estimate  of  the  orbit  is  ap¬ 
proaching  the  error  in  the  observations.  Under  this  assumption 


1/2 

improved 


where 


is  the  standard  deviation  of  the  observation. 


(A.  15) 


This  method  is  used  in  the  triple  filtering  mode,  updating  the  stand¬ 
ard  deviations  at  the  end  of  each  filtering.  It  appears  to  be  working 
quite  well. 


4 
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D.  REJECTION  CRITERION 


The  matrix  R  discussed  in  the  previous  section  is  composed  of  two 

T 

covariance  matrices.  Matrix  HPH  is  computed  by  the  program  itself. 
The  matrix  Q  is  an  input  item  for  the  first  filtering.  It  is  a  diagonal 
matrix  composed  of  the  variances  of  the  estimated  measurement  errors 
In  case  the  errors  are  not  known,  they  can  be  estimated  by  using  the 
triple  filtering  mode.  However,  an  initial  estimate  is  still  required. 

If  the  estimate  is  correct,  it  is  usually  standard  practice  to  reject  all 
data  that  exceed  3- sigma  accuracy.  However,  it  is  possible  that  the 
estimate  is  low,  in  which  case  legitimate  observations  may  be  rejected. 
To  avoid  this  possibility,  the  following  routine  is  used. 


At  each  observation  the  quantity  n 


y 


is  computed: 


1.25  I  AS  I 

,  xl/2 

(ry) 


(A.  16) 


where  AS  is  the  difference  between  the  measured  and  estimated  obser¬ 
vations,  and  r  is  the  corresponding  diagonal  element  in  the  R  matrix. 

The  n  's  are  stored  in  six  boxes  according  to  their  values 

N1  N2  N3  N4  N5  N6 
n  =  0  -  3  3-4  4-6  6-8  8-  11  11< 


At  the  end  of  each  filtering  phase,  the  quantities 
N.  N.  +  N9  N.  +...  +  N- 

0/  =  — ±-  Ui  =  — - - —  Ui  =  — - - - 

yl  N  *  2  N  '  •  •  •  N 


(A.  16a) 


are  computed,  where  Nj. . .  Ng  are  the  number  of  observations  falling 
into  a  particular  box.  N  is  the  total  number  of  observations. 

N  *  Nj  +  . . .  +  Ng  (A.  16b) 

The  ‘lys  are  then  tested  against  a  number  4*  to  determine  the  number 
of  sigmas  tolerated. 

4>f  <  4>  <.  4<i+1  (A.  16c) 

A  representative  value  of  4*  =  0.  95  was  established  by  experiment 
to  give  good  results.  The  number  of  sigmas  tolerated  in  the  next 
filtering  phase  is  found  from  the  following  table. 

i  =  1  2  3  4  5  6 

6  =  3  4  6  8  11  15  (A*  16d) 
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E.  REVOLUTION  NUMBER 


Revolution  number  is  defined  as  the  number  of  nodal  crossings  of 
the  equator  from  south  to  north.  An  initial  revolution  number  is  inputted 
in  the  program  and  the  first  revolution  is  added  at  the  first  nodal 
crossing.  To  save  computer  time,  the  nodal  crossing  is  determined 
from  a  second  order  nodal  period  and  finally  established  by  logical 
tests  which  are  described  in  the  revolution  number  subroutine. 

The  time  to  the  first  nodal  crossing  is  computed  from  the  input  or¬ 
bital  elements  as  follows. 


The  true  anomaly  of  the  node  is 

0-t  =  2r  -  oj 

N 


Then  the  eccentric  anomaly  is  obtained  from 


cos  8n  +  e 

Cos  EAT  =  g-—  ■■■"q 
N  1  +  e  cos  0 


-\ 


Sin  En  = 


(1  -  e2) 


1/2 


N 


sin0 


>  0  <En<  2» 


N 


1  +  e  cos  0 


N 


and  the  corresponding  mean  anomaly 
Mj^  =  Ejy  -  e  sin  Ej^ 

AMn  =  Mn  -  M 
where  0  <.  AM  S.  2w 


(A.  17) 

(A.  18a) 

(A. 18b) 

(A.  19) 
(A.  20) 


The  time  to  the  first  nodal  crossing  is 
AM 


tN  ' 


(?) 


TTT 


(A.  21) 


To  determine  the  revolution  number  at  any  arbitrary  time,  Tpr#  we 


obtain 


m0  = 


Tpr  "  Tin  ~  *N 
PN 


(A.  22) 
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where  TJn  is  the  initial  time  and  is  the  approximate  nodal  period 

PN=  IT  I1’!  J2  [S'4  -■in2i(4-|«2)]}  (A. 23) 

The  following  values  then  are  used  in  the  logical  tests  (see  sub¬ 
routine  PROUT). 


m  =  mQ  +  1 


|mj  =  absolute  value  of  fraction  of  m 
±m2  =  Integer  of  m 
±m^  =  integer  of 


(A.  23a) 
(A.  23b) 
(A.  23c) 
(A.  23d) 


F.  CONVERSION  FROM  POSITION  AND 
VELOCITY  VECTORS  TO  ORBITAL  ELEMENTS 


Given  are  the  components  of  the  position  and  velocity  vectors  x,  y, 
z,  x,  y,  z,  and  the  constants. 


Components  in  an  orbital  axes  system  (yw  =0,  zw  =  0)  are 
2  2 

=  (X2  +  y2  +  Z2) 

V2«  ^2-i2+y2+i2 
•  .  xx  +  yy  +  zz 


w 


w 


=  (v«  •  *1) 


1/2 


(A.  24) 
(A.  25) 
(A.  26) 

(A.  27) 


The  planar  orbital  elements  associated  with  the  equations  of  conic 
sections  in  Cartesian  form  (see  Appendix  of  Ref.  1)  are 

.2 


3 

>  j 

3 

j 

1 

II 

< 

(A- 28) 

B  =  x  x  y 

03  w  w 

(A-29) 

(x  y  )2 

r  03  03 

u  M- 

(A.  30) 
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and  eccentricity 


2  o  A/  * 
e  *  (A^  +  B  / 

Semimajor  axis  is  obtained  from 

-rS- 

1  -  e 


(A.  31) 


(A.  32) 


To  obtain  the  remaining  elements,  we  compute  the  direction  cosines  of 

the  x  axis 
(0 


5  s  JL.  „  =  JL  ?  a  z_ 
1  x, ,  '1  x,. #  S1  x,.» 


(A.  33) 


The  direction  cosines  of  the  zw  axis  are  obtained  by  taking  the  vector 
product  of  the  position  and  velocity  vectors 


Cg  =  y  z  -  z  y,  a  z  X  -  X  z,  ?g  !xy  -  yx 

(A.  34) 

d  -  («;*  +  ni*  +  p,»y/2 

(A.  35) 

ci 

?3  *  T>  V  d  '  r3  d 

(A.  36) 

Consequently,  the  inclination  is  given  by 

Bln  1  *  (l  -  | 

/  0  <  i  <  * 

cos  i  =  ? 3  ) 

(A.  37) 

The  right  ascension  of  the  ascending  node  is  obtained  from  the  vector 
product  of  the  normal  to  the  orbital  plane  (zw  axis)  and  the  normal  to 

the  reference  plane  (z  axis),  which  gives 


sin  O  =  s 


■  (■ 


?3+”3 


n  *  5„  * 


"  ('M) 


o  <  n  <  2t 


(A.  38) 
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The  angle  ^  Is  obtained  by  taking  first  the  scalar  product  of  the  xw 
axis  and  the  line  of  nodes,  and  then  the  xw  axis  and  an  axis  in  the 
orbital  plane  normal  to  the  line  of  nodes 


cos  +  *  ^  +  r?1  r)N 

sin  *  r  ?3  +  ?1  nN  *  ^3  W 

0  <  *l><  2* 

The  angle  0  is  obtained  from  the  relationship 


(A.  39) 


sin  0  «  -  §■ 
e 


cos  0  *  -  4 
e 


0  <  0  <  2» 


(A.  40) 


It  is  measured  from  the  x^  axis -positive  counterclockwise.  Then  the 
argument  of  perigee  is  simply 

w  «  +  +  0  0  <  u<  2*  (A.  41) 

The  eccentric  anomaly  is  obtained  from  the  familiar  relationships 
cos  0  +  e 


cos  E  * 


sin  E 


1  +  e  cos  0 


-  (1  -  a  ) 


1/2 


sin  0 


0  <  E0  <  2r 


(A.  42) 


1  +  e  cos  0 

From  which  ' 

M  *  E  -  e  sin  E 


(A.  43) 
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It  must  be  remembered  that  the  x  ,  y  .  z  axes  system  is  considered 

(i)  ll)  ll) 

here  as  an  inertial  axes  system  with  the  x  axis  pointing  at  the  satellite 
at  the  particular  instant  of  time.  u 

G.  CONVERSION  TO  POSITION  AND  VELOCITY  COORDINATES 

First,  we  find  the  transformation  equations  from  the  reference 
axes  system  x,  y,  z  to  the  orbital  axes  system  X,  Y,  Z.  The  reference 
axes  system  will,  normally,  be  the  equatorial  axes  system  with  x  axis 
towara  the  vernal  equinox,  z  axis  pointing  north  and  y  axis  in  the 
equatorial  plane  completing  a  right  hand  system. 

Starting  with  the  reference  system  we  rotate  about  z  axis  through 
the  angle  ft . 

Xj  *  cos  fix  +  sin  Gy  (A.  44a) 

yj  *  -  sin  Ox  +  cos  Gy  (A.  44b) 

Zj  *  z  (A.  44c) 

Next  the  system  is  rotated  about  the  new  x^  axis  through  the  angle  i 


x2,xl 

(A.  45a) 

y2  *  cos  i  yj  +  sin  i  Zj 

(A.  45b) 

z2  9  “  s*n  *  +  cos  i  zi 

(A.  45c) 

Finally,  we  rotate  the  system  about  the  z2  axis  through  the  angle  w 
X  *  cos  «  x2  +  sin  u  y2  (A.  46a) 

Y  *  “  sin  «  x2  +  cos  w  y2  (A.  46b) 

Z  *  ^2  (A.  46c) 

The  operation  can  be  written  in  matrix  form 


‘JT 

cos  <•>  sinw 

0 

1 

0  0 

cos  Q  sin  0  0 

X 

Y 

-  sin  <*>  cos  (*> 

0 

0 

cos  i  sin  i 

-  sin  ft  cosDO 

y 

Z 

•  • 

0  0 

1 

0 

-  sin  i  cos  i 

0  0  1 

z 

•  «■ 

(A.  47) 
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Performing  the  matrix  multiplication  we  obtain 

x 


Z 

Y 

z 


p  p  p 
x  y  z 

Qx  Qy  Qz 

Rx  Ry  Rz 


where 


Px  *  cos  <4  cos  O  -  sin  <4  sin  0  cos  i 


P  *  cos  (u  sin  Q  +  sin  <4  cos  ft  cos  i 

y 


P  *  sin  <4  sin  i 
z 


Qx  *  -  sin  w  cos  ft  -  cos  w  sin  ft  cos  i 


S 


*  -  sin  <4  sin  ft  +  cos  <4  cos  ft  cos  i 


Q  *  cos  <4  sin  i 
z 


R^  *  sin  ft  sin  i 


R  *  -  cos  ft  sin  i 

•7 


R  «  cos  i 
z 


Since  this  is  an  orthogonal  transformation 
x 


P  Q  R 

X  X  X 

P  Q  R 

y  y  y 

P„  Q  R 

Z  z  z 


X 

Y 

Z 


Next  we  transform  to  an  axes  system  where  the  new  x 
the  direction  of  the  satellite  at  time  t. 


Cl 


This  means  a  rotation  about  Z  axis  through  the  angle  9 
transformation  equations 


(A.  48) 


(A.  49) 

)  axis  is  in 
giving  the 


X 

1 

1 _ 

cos  0  sin  0  0 

XW 

y 

■ 

pyQyIV 

-  sin  0  cos  0  0 

y« 

z 

I 

C 

N 

N50 

1 _ 

0  0  1 

z“ 

Since 


(A.  50) 


A  *  -  e  cos  0 
B  *  -  e  sin  0 
C  *  a  (1  -  e2) 

we  obtain,  in  conjunction  with  Eqs  (A.  12),  (A.  13)  and  (A.  14)  of  Ref.  1, 


y« 


a  (1  -  e  ) 
1  +  e  cos  0 


0 

-  e  sin  0 

(1  +  e  cos 


i-^r2 

[a  (1  -  eZ)J 

6>  [a(l-e2)] 


1/2 


z  *  0 
w 


(A. 51a) 
(A.  51b) 

(A. 51c) 

(A. 51d) 
(A. 51e) 


Introducing  the  velocity  coordinates  and  multiplying  the  last  column 
vector  by  tne  matrix  in  Eqs  (A.  50),  we  obtain  a  column  matrix 


X 

a  (1  -  e2)  cos  0 

1  +  e  cos  0 

Y 

a  (1  -  e2)  sin  0 

1  +  e  cos  0 

Z 

= 

0 

• 

3T 

.  r_*_J 

La  (1  -  eZ)J 

• 

Y 

• 

(e  +  cos  0)  - ^ — k- 

La(l  -  e  jJ 

Z 

0 

(A.  52) 
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To  introduce  the  eccentric  anomaly,  we  use  the  relationships 


cos  6 


cos  E  -  e 
1  - e  cos  £ 


sin  8  « 


9  1/2 

-  (1  -  e *)  sin  E 

1  -  e  cos  E 


(A. 53a) 
(A. 53b) 


which  give 


f  a  (cos  E  -  e) 

”  5T“ 

,  1/2 

a  (1  -  e  )  sin  E 

Y 

0 

s 

0 

sin  E 

(!01/2 

• 

ST 

1  -  e  cos  E 

A 

cos  E_ 

0  (1 
1  -  e  cos  E  1 

-  e2)1'2  ( 

io1/2 

• 

Y 

0 

0 

(A-54) 


The  transformation  equations  in  matrix  form  can  now  be  written 
x  !  I  px  Qx  Rx  0  0  0 


PyQyRy0  0  0 

P2QZRZ°  0  0 

0  0  0  P  Q  R 

Jv  A  A 

0  0  0  PyQyRy 

o  o  o  PZQZRZ 


X 


X 

Y 

0 

• 

3T 

• 

Y 

0 

•  « 

(A.  55) 


which  can  be  simplified  to  the  following  form 
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(A.  56) 


The  advantage  of  forms  (A.  5b)  and  (A.  56)  is  that  the  elements  which 
determine  the  Keplerian  orbit  a,  e,  E  and  the  elements  which  determine 
the  orientation  of  the  orbit  ft,  i,  u  are  completely  separated. 


H.  COMPUTATION  OF  EPHEMERIS 

In  one  part  of  the  ephemeris  computations,  we  determine  the  satel¬ 
lite  longitude,  geocentric  and  geodetic  latitude  and  altitude  above  the 
reference  ellipsoid  from  the  position  and  velocity  vector  at  a  given 
time. 


where  Rg  is  the  Earth's  mean  equatorial  radius  and 

fE  *  (1  -  f)2  (A.  59) 

where  f  is  the  flattening. 

The  geodetic  latitude  is  obtained  by  a  successive  approximation 
technique  which  is  most  suitable  for  an  electronic  computer.  Com¬ 
puting  a  value 
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c  = 


2  ~\ 
1  +  tan  $ 

1  +fE  tan2  $  J 


(A.  60) 


we  substitute  tan  4>  for  tan  #  for  the  first  iteration.  Then  the  geodetic 
latitude  is 

h  +  C 
h  +fE  e 

where 


tan  4> 


(A.  61) 


♦  =  arctan 


An  improved  C  is  then  computed  using  tan  $  in  Eq  (A.  60)  and  the  process 
repeated  to  the  desired  accuracy.  Convergence  usually  occurs  in  one  or 
two  iterations. 

Satellite  longitude  is  determined  as  follows.  The  right  ascension  of 
the  satellite  is 

Xg  *  arctan  (jf^  (A.  62) 

The  right  ascension  of  the  Greenwich  meridian  is 

XG  *  1.  746647719  +  6.  30038809863056  d  +  0.  5064  x  10-14  d2  +  AX 

(A.  63) 

where  d  are  Julian  days  from  the  epoch  1950  January  1,  0h  UT  (see 
Ref.  1,  III-D).  The  satellite  longitude  is  then 

X  *  XG  “  xs  (A.  64) 
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